How to compute $I_n=\int_{-\infty}^{+\infty}\mathrm{d}x\frac{x^{2n}}{\cosh^2 x}$? I'd like to compute:
$$
I_n = \int_{-\infty}^{+\infty}\mathrm{d}x\frac{x^{2n}}{\cosh^2 x}.
$$
We have, quite easily:
$$
I_0 = \int_{-\infty}^{+\infty}\mathrm{d}x\frac{1}{\cosh^2 x}=\left[\tanh x\right]_{-\infty}^{+\infty}=2.
$$
So to begin with, I gave a try at
$$
I_1 = \int_{-\infty}^{+\infty}\mathrm{d}x\frac{x^{2}}{\cosh^2 x}.
$$
My idea was to use the following contour $\Gamma_{\varepsilon, M}$ to integrate on in the complex plane: the rectangle with vertices $M$, $M+i\frac{\pi}{2}$, $-M+i\frac{\pi}{2}$, $-M$, and an indent of radius $\varepsilon$ at $i\frac{\pi}{2}$, where $\cosh z=0$.
Then by Residues Theorem:
$$
\oint_{\Gamma_{\varepsilon, M}} \mathrm{d}z\frac{z^{2}}{\cosh^2 z} = 0.
$$
Problem is, when I make the calculations the integrand is still singular $O(\frac{1}{\varepsilon})$ on my semicircle. More explicitly:
$$
\int_{0}^{-\pi} \mathrm{d}\theta i\varepsilon e^{i\theta}\frac{(i\frac{\pi}{2}+\varepsilon e^{i\theta})^{2}}{\cosh^2 (i\frac{\pi}{2}+\varepsilon e^{i\theta})} = \int_{0}^{-\pi} \mathrm{d}\theta i\varepsilon e^{i\theta} \frac{-\frac{\pi^2}{4} + O(\varepsilon)}{\varepsilon^2e^{2i\theta}+O(\varepsilon^3)}.
$$
To avoid this problem I tried:
$$
\oint_{\Gamma_{\varepsilon, M}} \mathrm{d}z\frac{(z-i\frac{\pi}{2})^{2}}{\cosh^2 z}=0.
$$
But then I only get:
$$
\int \mathrm{d}x\frac{(x-i\frac{\pi}{2})^{2}}{\cosh^2x}-\int \mathrm{d}x\frac{x^{2}}{\cosh^2(x-i\frac{\pi}{2})}=0
$$
$$
\int \mathrm{d}x\left(\frac{x^2}{\cosh^2x}+\frac{x^2}{\sinh^2x}\right) = \frac{\pi^2}{2}
$$
since $\cosh(x+i\frac{\pi}{2}) = i\sinh x$.
 A: Note that:
$$\int _{-\infty }^{\infty }{\frac {{x}^{2\,n}}{ \cosh^2 \left( x \right) 
 }}\,{\mathrm{d}x}=2\,\int _{0}^{\infty }\!{\frac {{x}^{2\,n
}}{\cosh^2 \left( x \right) 
 }}\,{\mathrm{d}x}\tag{1}$$
then, consider the following integral for $\Re(s)>1$:
$$\begin{aligned}
\int _{0}^{\infty }{\frac {{x}^{s}}{ \cosh^2 \left( x \right) 
 }}{dx}&=-2\,\int _{0}^{\infty }\!{x}^{s}{\frac {\mathrm{d}}{\mathrm{d}x}}
 \left(   \frac{1}{ 1+{{\rm e}^{2\,x}}}\right)\, {\mathrm{d}x},\\
\mbox{integration by parts...}\quad&=2\,s\int _{0}^{\infty }\!{\frac {{x}^{s-1}}{1+{{\rm e}^{2\,x}}}}\,{\mathrm{d}x},\\
\mbox{partial fractions...}\quad&=2\,s\int _{0}^{\infty }\!2\,{\frac {{x}^{s-1}}{1-{{\rm e}^{4\,x}}}}\,{\mathrm{d}x}-2\,s\int _{0}^{\infty }\!{\frac {{x}^{s-1}}{1-{{\rm e}^{2\,x}}}}\,{\mathrm{d}x},\\
\mbox{rescale the variables...}\quad&=-s \left( {2}^{-2\,s+2}-{2}^{1-s} \right) \int _{0}^{\infty }\!{\frac {{x}^{s-1}}{-1+{{\rm e}^{x}}}}\,{\mathrm{d}x},\\
\mbox{geometric series...}\quad&=-s \left( {2}^{-2\,s+2}-{2}^{1-s} \right) \sum _{n=1}^{\infty } \left( \int _{0}^{\infty }\!{x}^{s-1}{{\rm e}^{-xn}}\,{\mathrm{d}x} \right),\\
\mbox{rescale the variable...}\quad&=-s \left( {2}^{-2\,s+2}-{2}^{1-s} \right) \left(\sum _{n=1}^{\infty } \frac{1}{n^s}\right)\left(\int _{0}^{\infty }\!{x}^{s-
1}{{\rm e}^{-x}}\,{\mathrm{d}x}\right),\\
\mbox{function definitions...}\quad&=-s \left( {2}^{-2\,s+2}-{2}^{1-s} \right) \Gamma \left( s \right) \zeta  \left( s \right),\\
\mbox{absorb the s into}\,\, \Gamma \mbox{...}\quad&=\left( -{2}^{-2\,s+2}+{2}^{1-s} \right) \zeta 
 \left( s \right) \Gamma  \left( s+1 \right) 
\end{aligned} \tag{2}$$
where $\zeta$ is the Riemann zeta function. Then, if $s=2n,\, n\in \mathbb{Z}$:
$$\Gamma  \left( 2\,n+1 \right) = \left( 2\,n \right) !,\quad\zeta 
 \left( 2\,n \right) ={\frac { \left( -1 \right) ^{n+1}B_{2n}  \left( 2\,\pi  \right) ^{2\,n}}{2\left( 2\,n \right) !}} \tag{3}$$
where $B$ denotes the Bernoulli number (or polynomial later) and thus:
$$\int _{-\infty}^{\infty }\!{\frac {{x}^{2\,n}}{ \cosh^2 \left( x \right) }}{dx}= \left( 2-2^{2-2n} \right)\left( -1 \right) ^{n+1}B_{2n}\, {\pi }^{2\,n} \tag{4}$$
or if you prefer:
$$\int _{-\infty}^{\infty }\!{\frac {{x}^{2\,n}}{ \cosh^2 \left( x \right) }}{dx}=-2B_{2n}\left(\frac{1}{2} \right) 
 \left( i\pi  \right) ^{2\,n} \tag{5}$$
A: Why not simply go up to $\Im{z} = \pi$ so that the rectangle $C$ has vertices $\pm R$ and $\pm R+i \pi$?  Then we consider
$$\oint_C dz \frac{z^{2 n+1}}{\cosh^2{z}}$$
which is equal to
$$\int_{-R}^R dx \frac{x^{2 n+1}}{\cosh^2{x}} + i \int_0^{\pi} dy \frac{(R+i y)^{2 n+1}}{\cosh^2{(R + i y)}} \\ + \int_R^{-R} dx \frac{(x+i \pi)^{2 n+1}}{\cosh^2{x}} + i \int_{\pi}^0 dy \frac{(-R+i y)^{2 n+1}}{\cosh^2{(-R + i y)}}$$
As $R \to \infty$, the second and fourth integrals vanish.  Thus the contour integral is equal to
$$\int_{-\infty}^{\infty} dx \frac{x^{2 n+1}-(x+i \pi)^{2 n+1}}{\cosh^2{x}}$$
Note that the highest power in the numerator is $x^{2 n}$, and that all odd powers vanish upon integration.  Thus the contour integral is expressible in terms of integrals of lower powers:
$$-i \pi (2 n+1)\int_{-\infty}^{\infty} dx \frac{x^{2 n}}{\cosh^2{x}} +i \sum_{k=0}^{n-1} (-1)^k \binom{2 n+1}{2 k} \pi^{2 (n-k)+1}\int_{-\infty}^{\infty} dx \frac{x^{2 k}}{\cosh^2{x}}$$
By the residue theorem, the contour integral is also equal to $i 2 \pi$ times the residue of the integrand at the pole $z=i \pi/2$, which is a double pole.  
$$\operatorname*{Res}_{z=i \pi/2} \frac{z^{2 n+1}}{\cosh^2{z}} = \lim_{z\to i \pi/2} \frac{d}{dz}\frac{(z-i \pi/2)^2 z^{2 n+1}}{\cosh^2{z}} = (-1)^{n+1} (2 n+1) \left ( \frac{\pi}{2}\right )^{2 n}$$
In order to verify this you can proceed as follows:
$$
\lim_{z\to i \pi/2} \frac{d}{dz}\frac{(z-i \pi/2)^2 z^{2 n+1}}{\cosh^2{z}} = \lim_{u\to 0} \frac{d}{dz}\frac{u^2 (u+i\pi/2)^{2 n+1}}{-\sinh^2{u}} = \lim_{u\to 0} (u+i\pi/2)^n\ \frac{-\left((2n+1)u^2+2(u+i\pi/2)u\right)2\sinh u+2(u+i\pi/2)u^2\cosh u}{\sinh^3u} = \lim_{u\to 0} (u+i\pi/2)^n\ \frac{-(2n+3)u^3-i\pi u +2u^3 + i\pi u + O(u^4)}{u^3+O(u^4)}= (-1)^{n+1}(2n+1)\left(\frac{\pi}{2}\right) ^{2n}.
$$
You then need to find the integrals of lower powers.  Actually, you can get around this by setting up a system of equations for each $k$ from $0$ to $n$.  Let
$$I_k = \int_{-\infty}^{\infty} dx \frac{x^{2 k}}{\cosh^2{x}}$$
and
$$R_k = (-1)^{k} (2 k+1) \left ( \frac{\pi}{2}\right )^{2 k}$$
Then we have
$$(2 n+1) I_n - \sum_{k=0}^{n-1} (-1)^k \binom{2 n+1}{2 k} \pi^{2 (n-k)} I_k = 2 R_n$$
$$(2 n-1) I_{n-1} - \sum_{k=0}^{n-2} (-1)^k \binom{2 n-1}{2 k} \pi^{2 (n-1-k)} I_k = 2 R_{n-1}$$
$$\cdots$$
$$3 I_1 - \pi^2 I_0 = 2 R_1$$
$$I_0 = 2 $$
A: $\newcommand{\+}{^{\dagger}}
 \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
 \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
 \newcommand{\dd}{{\rm d}}
 \newcommand{\down}{\downarrow}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
 \newcommand{\fermi}{\,{\rm f}}
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{{\rm i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\isdiv}{\,\left.\right\vert\,}
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\, #1 \,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
 \newcommand{\sech}{\,{\rm sech}}
 \newcommand{\sgn}{\,{\rm sgn}}
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
 \newcommand{\wt}[1]{\widetilde{#1}}$
$\ds{I_{n}\equiv\int_{-\infty}^{\infty}{x^{2n} \over \cosh^{2}\pars{x}}\,\dd x:\
     {\large ?}.\qquad n = 0,1,2,3,\ldots}$.

\begin{align}
\color{#c00000}{I_{n}}&=8\int_{0}^{\infty}{x^{2n}\expo{2x} \over \pars{\expo{2x} + 1}^{2}}\,\dd x
=-4\int_{x\ =\ 0}^{x\ \to\ \infty}x^{2n}\,\dd\pars{1 \over \expo{2x} + 1}
\\[3mm]&=2\delta_{n0} + 8n
\color{#00f}{\int_{0}^{\infty}{x^{2n - 1}\expo{-2x} \over 1 + \expo{-2x}}\,\dd x}
\end{align}

\begin{align}
&\color{#00f}{\int_{0}^{\infty}{x^{2n - 1}\expo{-2x} \over 1 + \expo{-2x}}\,\dd x}
=\sum_{k = 0}^{\infty}\pars{-1}^{k}\int_{0}^{\infty}
x^{2n - 1}\expo{-\pars{2k + 2}x}\,\dd x
\\[3mm]&=\sum_{k = 0}^{\infty}{\pars{-1}^{k} \over \pars{2k + 2}^{2n}}
\int_{0}^{\infty}x^{2n - 1}\expo{-x}\,\dd x
=-\,{\Gamma\pars{2n} \over 2^{2n}}
\sum_{k = 1}^{\infty}{\pars{-1}^{k} \over k^{2n}}
\\[3mm]&=-{\Gamma\pars{2n} \over 2^{2n}}
\bracks{{1 \over 2^{2n}}\sum_{k = 1}^{\infty}{1 \over k^{2n}}
-\sum_{k = 1}^{\infty}{1 \over \pars{2k - 1}^{2n}}}
\\[3mm]&=-{\Gamma\pars{2n} \over 2^{2n}}
\bracks{{1 \over 2^{2n}}\sum_{k = 1}^{\infty}{1 \over k^{2n}}
-\sum_{k = 1}^{\infty}{1 \over k^{2n}}
+ \sum_{k = 1}^{\infty}{1 \over \pars{2k}^{2n}}}
\\[3mm]&=-{\Gamma\pars{2n} \over 2^{2n}}\pars{{1 \over 2^{2n - 1}} - 1}
\sum_{k = 1}^{\infty}{1 \over k^{2n}}
={\Gamma\pars{2n} \over 2^{4n - 1}}\pars{2^{2n - 1} - 1}\zeta\pars{2n}
\end{align}

$$\color{#66f}{\large I_{n}}
\equiv\int_{-\infty}^{\infty}{x^{2n} \over \cosh^{2}\pars{x}}\,\dd x
=\color{#66f}{\large2\delta_{n0}
+
n{2^{2n - 1} - 1 \over 2^{4n - 4}}\,\Gamma\pars{2n} \zeta\pars{2n}}
$$

A: I would like to post an answer based on a slight variation of @Graham Hesketh's method, both because I find it a bit simpler and because it reflects the steps leading to the integrals temselves in the context where they appear.
For $t>0$ and $s>0$, we consider
$$
\int_0^\infty \frac{x^{s-1}}{e^{tx}-1}dx=\sum_{n=1}^\infty\int_0^\infty x^{s-1}e^{-ntx}dx=t^{-s}\sum_{n=1}^\infty n^{-s} \int_0^\infty x^{s-1}e^{-x}dx=\frac{\zeta(s)\Gamma(s)}{t^s}\,,
$$
and note that
$$
\int_0^\infty \frac{x^{s-1}}{e^{tx}-1}dx -\int_0^\infty\frac{x^{s-1}}{e^{tx}+1}dx=2\int_0^\infty \frac{x^{s-1}}{e^{2tx}-1}dx = 2^{1-s}\int_0^\infty\frac{x^{s-1}}{e^{tx}-1}dx
$$
so, by comparison,
$$
\int_0^\infty \frac{x^{s-1}}{e^{tx}+1}dx=\frac{1-2^{1-s}}{t^s}\zeta(s)\Gamma(s)\,.
$$
The derivative of both sides of this equation, with respect to $t$, yields
$$
\int_0^\infty \frac{x^s}{\cosh^2(\frac{tx}{2})}dx = \frac{4s(1-2^{1-s})}{t^{s+1}}\zeta(s)\Gamma(s)\,.
$$ 
In particular, for $s=2n$ with $n\in\mathbb N$, and $t=2$ we retrieve the integrals we wanted
$$
I_{n}=\int_{-\infty}^{+\infty} \frac{x^{2n}}{\cosh^2 x}dx = 2^{1-2n}(1-2^{1-2n})\,4n\,\zeta(2n)\Gamma(2n)\,.
$$
E.g.
$$
I_0=2\,,\qquad I_1=\frac{\pi^2}{6}\,,\qquad I_2=\frac{7\pi^4}{120}\,.
$$
(The case $n=0$ requires formally an analytic continuation $\epsilon \Gamma(\epsilon)=1+\mathcal O(\epsilon)$ and $\zeta(0)=-1/2$).
