We first notice that
\begin{align*}
\int_{0}^{\infty}\frac{\operatorname{gd}(x)}{e^x-1}\,dx
&= \int_{0}^{\infty} \frac{1}{e^x-1} \left( \int_{0}^{x} \frac{dy}{\cosh y} \right) \, dx \\
&= \int_{0}^{\infty} \frac{1}{\cosh y} \left( \int_{y}^{\infty} \frac{dx}{e^x - 1} \right) \,d y \\
&= -2 \int_{0}^{\infty} \frac{\log(1 - e^{-y})}{e^y + e^{-y}} \, dy \\
&=-2\int_{0}^{\frac{\pi}{4}}\log(1-\tan\theta)\,d\theta \tag{$e^{-x}=\tan\theta$}.
\end{align*}
The last integral is our starting point. We introduce two tricks to evaluate this.
Step 1. Notice that $\tan(\frac{\pi}{4}-\theta)=\frac{1-\tan\theta}{1+\tan\theta}$. So by the substitution $\theta \mapsto \frac{\pi}{4}-\theta$, it follows that
$$ \int_{0}^{\frac{\pi}{4}}\log(1+\tan\theta)\,d\theta
= \int_{0}^{\frac{\pi}{4}}\log\left(\frac{2}{1+\tan\theta}\right)\,d\theta $$
and hence both integrals have the common value $\frac{\pi}{8}\log 2$. Applying the same idea to our integral, it then follows that
\begin{align*}
-2\int_{0}^{\frac{\pi}{4}}\log(1-\tan\theta)\,d\theta
&= -2\int_{0}^{\frac{\pi}{4}}\log\left(\frac{2\tan\theta}{1+\tan\theta}\right)\,d\theta \\
&= -2\int_{0}^{\frac{\pi}{4}}\log\tan\theta \, d\theta - \frac{\pi}{4}\log 2.
\end{align*}
Step 2. In order to compute the last integral, we notice that for $\theta\in\mathbb{R}$ with $\cos\theta\neq0$, we have
\begin{align*}
-\log\left|\tan\theta\right|
&= \log\left|\frac{1+e^{2i\theta}}{1-e^{2i\theta}}\right|
= \operatorname{Re} \log\left(\frac{1+e^{2i\theta}}{1-e^{2i\theta}}\right) \\
&= \operatorname{Re}\left( \sum_{n=1}^{\infty} \frac{1+(-1)^n}{n} e^{2in\theta} \right) \\
&= \sum_{k=0}^{\infty} \frac{2}{2k+1}\cos(4k+2)\theta.
\end{align*}
So by term-wise integration, we obtain
\begin{align*}
-2\int_{0}^{\frac{\pi}{4}}\log\tan\theta \, d\theta
&= \sum_{k=0}^{\infty} \frac{4}{2k+1} \int_{0}^{\frac{\pi}{4}} \cos(4k+2)\theta \, d\theta \\
&= 2 \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)^2}
= 2K,
\end{align*}
where $K$ is the Catalan's constant.
\begin{align}\mathcal{I}&=\frac{1}{2}\int_0^\infty \frac{x^n}{\cosh(x)}\mathrm dx\\&=\int_0^\infty \frac{x^n}{e^x+e^{-x}}\mathrm dx\\&=\int_0^\infty \frac{x^ne^{-x}}{1+e^{-2x}}\mathrm dx.\end{align}
Proposition: $$\int_0^\infty \frac{x^{s-1}e^{-x}}{1+e^{-2x}}\mathrm dx=\beta(s)\Gamma(s),$$ where $\beta(s)$ is the Dirichlet beta function defined as $\beta(s)=\sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)^s}$.
Proof:\begin{align}\int_0^\infty x^{s-1}\left(\frac{e^{-x}}{1-(-e^{-2x})}\right)\mathrm dx&=\int_0^\infty x^{s-1}\left(\sum_{k=0}^\infty (-1)^ke^{-(2k+1)x}\right)\mathrm dx \text{, using geometric series}\\&=\sum_{k=0}^\infty (-1)^k\int_0^\infty x^{s-1} e^{-(2k+1)x}\mathrm dx\\&=\sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)^s}\displaystyle\int_0^\infty x^{s-1}e^{-x}\mathrm dx\text{, substituting $(2k+1)x\mapsto x$}\\&=\beta(s)\Gamma(s). \end{align}
Therefore,
$$\mathcal{I}=\beta(n+1)\Gamma(n+1).$$
We can evaluate for different values of $n\ge 0$.
For $n=1$, $\frac{1}{2}\int_0^\infty x\operatorname{sech}(x)=\beta(2)\Gamma(2)=G$, where $G$ is Catalan's constant.
For $n=2$, $\int_0^\infty x^2\operatorname{sech}(x)=2\beta(3)\Gamma(3)=\frac{\pi^3}{8}$
Best Answer
Integrate as follows
\begin{align} &\int_0^{\infty}\frac{\operatorname{sech}(\pi x)}{1+4x^2}\>dx\\ \overset{t=2\pi x }=& \frac\pi2 \int_{-\infty}^{\infty}\frac{e^{\frac t2}}{(e^t+1)(\pi^2+t^2)}dt\\ =& \frac\pi2 \int_{-\infty}^{\infty}\frac{e^{\frac t2}}{e^t+1}Re\left(\frac1\pi \int_0^\infty e^{-(\pi-i t)y }dy\right)dt\\ =& \frac12 Re\int_{0}^{\infty}e^{-\pi y} \left(\int_{-\infty}^\infty \frac{e^{a t}}{e^t+1}dt \right)dy \>\>\>\>\>\>\>a= iy+\frac12\\ \overset{x=e^t}=& \frac12 Re\int_{0}^{\infty}e^{-\pi y} \left(\int_{0}^\infty \frac{x^{a-1}}{x+1}dx \right)dy \\ =& \frac12Re \int_{0}^{\infty}e^{-\pi y}\pi\csc(\pi a)\,dy \\ =& \int_{0}^{\infty} \frac\pi{e^{2\pi y}+1}dy \overset{t=e^{-2\pi y}}=\frac12\int_0^1 \frac1{1+t}dt\\ =&\frac12\ln2 \end{align}