We can use contour integration, but we have to be a little careful with it.

Begin by noting that the integrand is an even function, so we shall evaluate

$I=\int_{-\infty}^{+\infty}\dfrac{\sinh x}{x\cosh^2x}~dx$

and take half of this result. in this way the contour accesses infinity at both ends of the desired integration range, which facilitates using contour integration.

A semicircular contour that is often used with this doubly infinite integration range is problematic because the integrand will not be $o(1/z)$ over the entire semicircle. We violate the small-$o$ relation at the imaginary axis where there are infinitely many singularities.

We therefore design a more refined contour: a rectangle with corners at $-N, +N, +N+(2\pi N)i, -N+(2\pi N)i$. Traversing the rectangle counterclockwise and applying the Residue Theorem then gives

$I(1) + I(2) + I(3) + I(4)=2\pi i\sum _{k=0}^{2N} R[(2k+1)\pi i/2]$

where the integration path for $I(1)$ runs from $-N$ to $+N$, the path for $I(2)$ runs from $+N$ to $+N+(2\pi N)i$, and so on around the remaining sides of the rectangle for $I(3)$ and $I(4)$. Then the required integral is the limit of $I(1)$ as $N\to\infty$, which matches with the sum of residues provided that $I(2),I(3),I(4)$ go to zero in this limit.

To verify the zero limit for $I(2)$ render (limits represent the asymptotic behavior as $N\to\infty$):

$|\sinh(N+iy)|=\sqrt{\sinh^2N+\sin^2y}\to e^{N/2}$

$|\cosh(N+iy)|=\sqrt{\cosh^2N-\sin^2y}\to e^{N/2}$

Thus since the length of the side is $2\pi N$, the Triangle Inequality guarantees

$|I(2)|\le\dfrac{2\pi\sqrt{\sinh^2N+\sin^2y}}{(\cosh^2N-\sin^2y)}\to 2\pi e^{-N/2}\to0$

A similar analysis eliminates $I(4)$ as $N\to\infty$.

For $I(3)$ we use the fact that the hyperbolic functions are periodic, thus for the selected values of $y$ they match the real-argument values at $x$. Thereby

$|\dfrac{\sinh z}{z\cosh^2 z}|=|\dfrac{\sinh x}{z\cosh^2 x}|<\dfrac{1}{y\cosh x}<[1/(2\pi N)]e^{-|x|/2}$

$|I(3)|<[1/(2\pi N)]\int_{-N}^Ne^{-|x|/2}dx\to0$ (the integral is bounded as $N\to\infty$.)

So we have

$I=\lim_{N\to\infty}(I(1))=2\pi i\sum _{k=0}^\infty R[(2k+1)\pi i/2]$

and it remains to find the residues.

To find these residues $R[(2k+1)\pi i/2]$, it is convenient to define a difference parameter $\delta$ at each of these singularities. Since the singularities are second-order poles, we need to carry the Laurent series expansion to two terms.

Therefore, at each singularity $z=(2k+1)\pi i/2$, we render the following Laurent series:

$\sinh z = (-1)^ki+0\delta+O(\delta^2)$

$\dfrac{1}{\cosh^2 z} = 1/\delta^2+0/\delta+O(1)$

$(1/z)=\dfrac{-2i}{(2k+1)\pi}-\dfrac{4}{(2k+1)^2\pi^2}\delta+O(\delta^2)$

And upon multiplying theses:

$\dfrac{\sinh z}{z\cosh^2z}=\dfrac{2(-1)^k(2k+1)}{\pi\delta^2}\color{blue}{-\dfrac{4(-1)^ki}{(2k+1)^2\pi^2\delta}}+O(1)$

from which we read off the residues

$R[(2k+1)\pi i/2]=\dfrac{-4(-1)^ki}{(2k+1)^2\pi^2}$

So

$I=\int_{-\infty}^{+\infty}\dfrac{\sinh x}{x\cosh^2x}~dx=(8/\pi)\color{blue}{\sum _{k=0}^\infty \dfrac{(-1)^k}{(2k+1)^2}}$

where the blue series directly defines the Catalan constant $G$. Thus $I=8G/\pi$ and the halved value in the original problem becomes

$\int_0^\infty\dfrac{\sinh x}{x\cosh^2x}~dx=4G/\pi.$