Find the integral of $\frac{x^n}{(\exp(x)+1)(\exp(-x)+1)}$ I need to find the integral $$I = \int_{-\infty}^\infty \frac{x^n}{(\exp(x)+1)(\exp(-x)+1)}dx.$$
I think this should be possible with the residue theorem, but I got stuck.
So first of all the integral is zero if n is odd, since then it is an odd integrand integrated over a symmetric interval.
My choice for the complex contour is a rectangle with sidelength $2R$ (on the real axis) and the other $2\pi i$, because the denominator is the same for $x \rightarrow x+2\pi i$ and $I = \lim_{R\rightarrow \infty} \int_{\gamma_1}f(z)dz$.

*

*I define the contour $\gamma = \gamma_1 + \gamma_2 + \gamma_3 + \gamma_4$, where $\gamma_1$ and $\gamma_3$ are the two long, horizontal sides (with length $2R$) and $\gamma_2$ and $\gamma_4$ the two vertical ones. Then there is one pole of order 2 inside the region of integration at $z_0 = i\pi$.


*$\int_{\gamma_2} f(z)dz + \int_{\gamma_4}f(z)dz = \int_0^{2\pi i} \frac{(t+R)^n}{(\exp(t+R)+1)(\exp(-t-R)+1)}dt - \int_0^{2\pi i} \frac{(-t+2\pi i-R)^n}{(\exp(-t-R)+1)(\exp(+t+R)+1)}dt = 0 + 0$ individually when $R\rightarrow \infty$.


*Now this is where I got stuck. I want to relate the integral over $\gamma_3$ to the integral over $\gamma_1$. For this I tried to rearange
$$\begin{align}
\int_{\gamma_3}f(z)dz  &= -\int_{-R}^{R}\frac{(-t+2\pi i)^n}{(\exp(t) + 1)(\exp(-t)+1)}dt \\
\\
 &= -\sum_{k=0}^n \binom{n}{k}\cdot\int_{-R}^R\frac{(-t)^k(2\pi i)^{n-k}}{(\exp(t)+1)(\exp(-t)+1)}dt \\
\\
&=  
-\sum_{k=0}^n \binom{n}{k}\cdot(2\pi i)^{n-k}\cdot I_k
\end{align}$$
where I defined $I_k = \int_{\gamma_1}f(z)dz$.
Now I don't know how to proceed from here. Is there a way I can solve this recursive formula or do I is there maybe a completely different way of solving this?
 A: $$J(n)=\int_{-\infty}^\infty \frac{x^n}{(e^x+1)(e^{-x}+1)}dx=\frac{(2\pi)^n\pi}{2}\int_{-\infty}^\infty\frac{x^n}{\cosh^2\pi x}dx$$
$$=\frac{(2\pi)^n\pi}{2}\frac{\partial^n}{\partial \alpha^n}\Big|_{\alpha=0}\int_{-\infty}^\infty\frac{e^{\alpha x}}{\cosh^2\pi x}dx=\frac{(2\pi)^n\pi}{2}\frac{\partial^n}{\partial \alpha^n}I(\alpha)\Big|_{\alpha=0}$$
The integral $\displaystyle I(\alpha)=\int_{-\infty}^\infty\frac{e^{\alpha x}}{\cosh^2\pi x}dx$ is well-known and can be evaluated, for example, via complex integration.
Consider the following rectangular contour

and the integral along it gives
$$\oint\frac{e^{\alpha z}}{\cosh^2\pi z}dz=I(\alpha)+[1]-I(\alpha)e^{i\alpha}+[2]=2\pi i\underset{z=\frac{i}{2}}{\operatorname{Res}}\frac{e^{\alpha z}}{\cosh^2\pi z}$$
because we have a single second order pole inside the contour. Integrals $[1]$ and $[2]$ tend to zero, as $R\to\infty$.
Near $z=\frac{i}{2}$ we have
$$\frac{e^{\alpha \frac{i}{2}+\epsilon}}{\cosh^2\pi(\frac{i}{2}+\epsilon)}=\frac{e^{i\frac{\alpha}{2}}}{-\pi^2\epsilon^2}\big(1+\alpha\epsilon+O(\epsilon^2)\big)$$
so, the residue is $\displaystyle -\,\frac{\alpha e^{i\frac{\alpha}{2}}}{\pi^2}$
$$I(\alpha)\big(1-e^{i\alpha}\big)=-\frac{2i\alpha}{\pi}e^{i\frac{\alpha}{2}}\,\,\Rightarrow\,\,\boxed{\,I(\alpha)=\frac{\alpha}{\pi\sin\frac{\alpha}{2}}\,}$$
$$J(n)=\frac{(2\pi)^n\pi}{2}\frac{\partial^n}{\partial \alpha^n}I(\alpha)\Big|_{\alpha=0}=\boxed{\,\pi^n\frac{d^n}{dx^n}\bigg(\frac{x}{\sin x}\bigg)\Big|_{x=0}\,}$$
Following the suggestion by @Gary (https://dlmf.nist.gov/4.19#E4), the answer can be expressed via Bernoulli numbers (probably, this is the most compact form):
$$\boxed{\,\,I(2k)= (-1)^{k-1}2\pi^{2k}(2^{2k-1}-1)B_{2k},\,\,\text{and}\,\,I(2k+1)=0;\,\, k=0, 1, 2, ...\,\,}$$

Quick check: at $n=0\,\, J(0)=\pi^0=1$
On the other hand, $\displaystyle J(n=0)=\int_0^\infty\frac{dt}{\cosh^2t}=\tanh t\Big|_0^\infty=1$
A: At first we see that
\begin{gather*}
 I(n)=\int_{-\infty}^{\infty}\dfrac{x^n}{(e^x+1)(e^{-x}+1)}\,\mathrm{d}x= \int_{-\infty}^{\infty}\dfrac{x^ne^x}{(e^x+1)^2}\,\mathrm{d}x =\\[2ex]
 \int_{0}^{\infty}\dfrac{x^ne^x}{(e^x+1)^2}\,\mathrm{d}x+\int_{0}^{\infty}\dfrac{(-x)^ne^{-x}}{(e^{-x}+1)^2}\,\mathrm{d}x=\\[2ex]
 (1+(-1)^n)\int_{0}^{\infty}\dfrac{x^ne^x}{(e^x+1)^2}\,\mathrm{d}x.
\end{gather*}
Apparently
\begin{equation*}
 I(n)=0 
\end{equation*}
if $n$ is an odd and positive integer.
If $k>0$ then integration by parts yields
\begin{equation*}
 I(2k)=-2\underbrace{\left[\dfrac{1}{e^x+1}x^{2k}\right]_{0}^{\infty}}_{= 0}+2\cdot2k\int_{0}^{\infty}\dfrac{x^{2k-1}}{e^x+1}\,\mathrm{d}x=4k\int_{0}^{\infty}\dfrac{x^{2k-1}}{e^x+1}\,\mathrm{d}x. \tag{1}
\end{equation*}
However,
\begin{equation*}
 \int_{0}^{\infty}\dfrac{x^{s-1}}{e^x-1}\,\mathrm{d}x =\Gamma(s)\zeta(s).
\end{equation*}
See
https://en.wikipedia.org/wiki/Riemann_zeta_function
We get that
\begin{gather*}
 \int_{0}^{\infty}\dfrac{x^{s-1}}{e^x+1}\,\mathrm{d}x = \int_{0}^{\infty}x^{s-1}\left(\dfrac{1}{e^x+1}-\dfrac{1}{e^x-1}\right)\,\mathrm{d}x+\int_{0}^{\infty}\dfrac{x^{s-1}}{e^x-1}\,\mathrm{d}x =\\[2ex]
 -2\int_{0}^{\infty}\dfrac{x^{s-1}}{e^{2x}-1}\,\mathrm{d}x+\int_{0}^{\infty}\dfrac{x^{s-1}}{e^x-1}\,\mathrm{d}x=\\[2ex]
 (-2^{1-s}+1)\int_{0}^{\infty}\dfrac{x^{s-1}}{e^x-1}\,\mathrm{d}x=(1-2^{1-s})\Gamma(s)\zeta(s).
\end{gather*}
If we use that formula in $(1)$ we get
\begin{equation*}
 I(2k)=4k(1-2^{1-2k})\Gamma(2k)\zeta(2k)=2(1-2^{1-2k})(2k)!\zeta(2k)\tag{2}.
\end{equation*}
Since $\zeta(0)=-\dfrac{1}{2}$ formula $(2)$ is true also for $k=0.$
If we want to express the value of the integral using Bernoulli numbers we can use the formula
\begin{equation*}
 \zeta(2k)=\dfrac{(-1)^{k+1}B_{2k}(2\pi)^{2k}}{2(2k)!}
\end{equation*}
where $B_{2k}$ is the $2k$-th Bernoulli number.
See again
https://en.wikipedia.org/wiki/Riemann_zeta_function
Via $(2)$ we finally get
\begin{equation*}
 I(2k)=(-1)^{k-1}2\pi^{2k}(2^{2k-1}-1)B_{2k}.
\end{equation*}
A: 
$$J(n)=\int_{-\infty}^\infty \frac{x^n}{(e^x+1)(e^{-x}+1)}dx=\frac{(2\pi)^n\pi}{2}\int_{-\infty}^\infty\frac{x^n}{\cosh^2\pi x}dx$$
$$=\frac{(2\pi)^n\pi}{2}\frac{\partial^n}{\partial \alpha^n}\Big|_{\alpha=0}\int_{-\infty}^\infty\frac{e^{\alpha x}}{\cosh^2\pi x}dx=\frac{(2\pi)^n\pi}{2}\frac{\partial^n}{\partial \alpha^n}I(\alpha)\Big|_{\alpha=0}$$

Based on the work from @Svyatoslav , an alternative method to evaluate $I(\alpha)$,
$$I(\alpha)=\int_{-\infty}^\infty\frac{e^{\alpha x}}{\cosh^2\pi x}dx=\int_{-\infty}^\infty\frac{e^{\alpha x}}{\left( \frac{e^{\pi x}+e^{-\pi x}}{2} \right)^2}dx=\int_{-\infty}^\infty\frac{4e^{(\alpha+2\pi)x}}{(e^{2\pi x}+1)^2}dx$$
Let: $t=e^{2\pi x},$
$$I(\alpha)=\frac{2}{\pi}\int_{0}^\infty\frac{t^{\frac{\alpha}{2\pi}}}{(1+t)^2}dt$$
Now, we can use beta function:
$$\text{B}(x,y)=\int_0^\infty \frac{t^{x-1}}{(1+t)^{x+y}}dt$$
So we get:
$$I(\alpha)=\frac{2}{\pi}\cdot\text{B}(1+\frac{\alpha}{2\pi},1-\frac{\alpha}{2\pi})=\frac{2}{\pi}\cdot\frac{\Gamma(1+\frac{\alpha}{2\pi})\cdot\Gamma(1-\frac{\alpha}{2\pi})}{\Gamma(2)}$$
Note:
$$\Gamma(2)=1,~ \text{and}~~~\Gamma(1+\frac{\alpha}{2\pi})=\frac{\alpha}{2\pi}\cdot\Gamma(\frac{\alpha}{2\pi})$$
So we get:
$$I(\alpha)=\frac{2}{\pi}\cdot\frac{\alpha}{2\pi}\cdot\Gamma(\frac{\alpha}{2\pi})\cdot\Gamma(1-\frac{\alpha}{2\pi})$$
Finally, use Euler's reflection formula:
$$\Gamma(\frac{\alpha}{2\pi})\cdot\Gamma(1-\frac{\alpha}{2\pi})=\frac{\pi}{\sin(\pi\cdot\frac{\alpha}{2\pi})}$$
So we are done for $I(\alpha)$:
$$I(\alpha)=\frac{\alpha}{\pi\cdot\sin(\frac{\alpha}{2})}$$
