How to compute $I_n=\int_{-\infty}^{+\infty}\mathrm{d}x\frac{x^{2n}}{\cosh^2 x}$?

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 $$

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}$$

$\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}} $$