Use the exponential definition, $\operatorname{sech}(x)=\frac{2}{e^{x}+e^{-x}}$, to rearrange the integrand.
$$\begin{aligned}\operatorname{sech}^2(x-a)\operatorname{sech}^2(x+a)&=\left(\frac{2}{e^{x-a}+e^{-(x-a)}}\right)^2\left(\frac{2}{e^{x+a}+e^{-(x+a)}}\right)^2
\\
&=\frac{16}{\left(e^{2x}+e^{-2x}+c\right)^2}
\end{aligned}$$
where $c=e^{2a}+e^{-2a}=\frac{2}{\operatorname{sech}(2a)}$. The integrand is an even function, so $f(a)=32\int_{0}^{+\infty}\frac{1}{\left(e^{2x}+e^{-2x}+c\right)^2}\ \mathrm{d}x$. Substitute $t=e^{2x}$, $\mathrm{d}x=\frac1{2t}\mathrm{d}t$.
$$\begin{aligned}f(a)&=32\int_1^{+\infty}\frac{1}{(t+1/t+c)^2}\cdot\frac{1}{2t}\ \mathrm{d}t
\\
&=16\int_1^{+\infty}\frac{t}{(t^2+ct+1)^2}\ \mathrm{d}t
\end{aligned}$$
It then suffices to use the standard integrals of $\frac{x}{(Ax^2+Bx+C)^n}$ and $\frac{1}{Ax^2+Bx+C}$. To avoid confusion (since unfortunately we've chosen the same letters), I'll only refer to the constants from the integral table in capitals.
$$\begin{aligned}f(a)&=16\left[-\frac{ct+2}{(4-c^2)(t^2+ct+1)}-\frac{c}{4-c^2}\int\frac{1}{(t^2+ct+1)}\ \mathrm{d}t\right]_1^{+\infty}
\\
&=16\left[-\frac{ct+2}{(4-c^2)(t^2+ct+1)}-\frac{c}{4-c^2}\left(\frac{1}{\sqrt{c^2-4}}\ln\left|\frac{2t+c-\sqrt{c^2-4}}{2t+c+\sqrt{c^2-4}}\right|\right)\right]_1^{+\infty}
\end{aligned}
$$
Note that $4-c^2<0$ except when $a=0$. Hence, we used the integral in the table for the case $4AC-B^2<0$.
As $t\to+\infty$, the $t^2$ in the denominator of the first term dominates, so the term tends to $0$. Likewise, the argument of the logarithm tends to $1$, so the second term tends to $0$. Hence, we just substitute $t=1$ and remember to take the negative sign since it's the integral's lower bound.
$$
\begin{aligned}
f(a)&=-16\left[-\frac{c+2}{(4-c^2)(1+c+1)}-\frac{c}{4-c^2}\left(\frac{1}{\sqrt{c^2-4}}\ln\left|\frac{2+c-\sqrt{c^2-4}}{2+c+\sqrt{c^2-4}}\right|\right)\right]
\\
&=
\frac{16}{(4-c^2)}\left(1+\frac{c}{\sqrt{c^2-4}}\ln\left(\frac{c-\sqrt{c^2-4}}{2}\right)\right)
\\
\end{aligned}$$
Note that $c>0$, which allowed us to simplify the absolute value in the logarithm.
Now, to find $f(a)$ in terms of $a$, we substitute back $c=e^{2a}+e^{-2a}$ and find
$$
\begin{aligned}
f(a)&=\frac{16}{-(e^{4a}-2+e^{-4a})}\left(1+\frac{e^{2a}+e^{-2a}}{\sqrt{(e^{4a}-2+e^{-4a})}}\ln\left(\frac{e^{2a}+e^{-2a}-\sqrt{(e^{4a}-2+e^{-4a})}}{2}\right)\right)
\\
&=\frac{-16}{(e^{2a}-e^{-2a})^2}\left(1+\frac{e^{2a}+e^{-2a}}{e^{2a}-e^{-2a}}\ln\left(\frac{e^{2a}+e^{-2a}-e^{2a}+e^{-2a}}{2}\right)\right)
\\
&=\frac{16}{(e^{2a}-e^{-2a})^2}\left(2a\frac{e^{2a}+e^{-2a}}{e^{2a}-e^{-2a}}-1\right)
\\
&=4\left(2a\frac{e^{2a}+e^{-2a}}{e^{2a}-e^{-2a}}-1\right)\left(\frac{2}{e^{2a}-e^{-2a}}\right)^2
\\
&=4\left(2a\coth(2a)-1\right)\operatorname{csch}^2(2a)
\end{aligned}$$
Where, in the second line, we have used $(e^{4a}-2+e^{-4a})=(e^{2a}-e^{-2a})^2$. Hence, the solution to the integral agrees with the one given by Wolfram Alpha.
Here is some similarity between your attempt and a complex approach. Substitute $z=dx$ and write $$I=\frac d2\Re\int_\Bbb R\frac{\exp(-A/(z^2+B))}{z^2+C}e^{iz}\,dz$$ where $A=ad^2$, $B=bd^2$ and $C=cd^2$. Take a semi-circular contour on the upper-half plane so that the line integral on $[-R,R]$ gives us $I$. The integral along the arc $Re^{it}$ $(R\to\infty;0\le t\le\pi)$ tends to zero by Jordan's lemma, so $$I=\frac d2\cdot2\pi\Re i(\operatorname{Res}(f,i\sqrt C)+\operatorname{Res}(f,i\sqrt B))$$ by Cauchy's theorem ($f$ denotes the integrand of $I$). We have \begin{align}\operatorname{Res}(f,i\sqrt C)&=\lim_{z\to i\sqrt C}\frac{e^{-A/(z^2+B)}}{z+i\sqrt C}e^{iz}=\frac{e^{A/(C-B)}}{2i\sqrt C}e^{-\sqrt C}=\frac{e^{a/(c-b)}}{2id\sqrt c}e^{-d\sqrt c}\\\operatorname{Res}(f,i\sqrt B)&=\operatorname{Res}\left(\frac{e^{i(z+i\sqrt B)-A/((z+i\sqrt B)^2+B)}}{(z+i\sqrt B)^2+C},0\right)=e^{-d\sqrt b}\operatorname{Res}\left(\frac{e^{iz}e^{-A/(z(z+2id\sqrt b))}}{(z+id\sqrt b)^2+cd^2},0\right)\end{align} so (note that as expected, the first term is the same as what you have computed) $$I=\frac{\pi e^{a/(c-b)-d\sqrt c}}{2\sqrt c}-\pi de^{-d\sqrt b}\Im\operatorname{Res}\left(\frac{e^{iz}e^{-A/(z(z+2id\sqrt b))}}{(z+id\sqrt b)^2+cd^2},0\right).$$
Best Answer
By http://dlmf.nist.gov/10.32.E8, we have \begin{align*} \int_0^{ + \infty } {e^{ - bx} \sqrt {x(x + a)} dx} & \mathop = \limits^{x = at} a^2 \int_0^{ + \infty } {e^{ - abt} \sqrt {t(t + 1)} dt} \\ & \! \mathop = \limits^{t = \frac{{w - 1}}{2}} \frac{{a^2 }}{4}e^{\frac{{ab}}{2}} \int_1^{ + \infty } {e^{ - \frac{{ab}}{2}w} \sqrt {w^2 - 1} d} w =\frac{a}{{2b}}e^{\frac{{ab}}{2}} K_1 \!\left( {\frac{{ab}}{2}} \right). \end{align*} This result is valid whenever $\Re a>0$ and $\Re b>0$.