Computing $\int_0^\infty\mathrm{d} x\frac{x}{e^x+1}$ with contour integration Let's set:
$$
\int_0^\infty\mathrm{d}x\frac{x}{e^x+1}=I.
$$
I would like to compute it using, presumably, the methods of complex analysis and contour integration.
 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{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}
 \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]{\,#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}$
$\large\tt\mbox{Just another one !!!}$

In the first step we integrate by parts:
  \begin{align}
&\color{#00f}{\large\int_{0}^{\infty}{x \over \expo{x} + 1}\,\dd x}
=\int_{0}^{\infty}\ln\pars{1 + \expo{-x}}\,\dd x=
\int_{0}^{\infty}\sum_{\ell = 1}^{\infty}{\pars{-1}^{\ell + 1} \over \ell}\,
\expo{-\ell x}\,\dd x
\\[3mm]&=\sum_{\ell = 1}^{\infty}{\pars{-1}^{\ell + 1} \over \ell}\
\overbrace{\int_{0}^{\infty}\expo{-\ell x}\,\dd x}^{\ds{1 \over \ell}}
=\sum_{\ell = 1}^{\infty}{\pars{-1}^{\ell + 1} \over \ell^{2}}
=\sum_{\ell = 0}^{\infty}
\bracks{{1 \over \pars{2\ell + 1}^{2}} - {1 \over \pars{2\ell + 2}^{2}}}
\\[3mm]&=\bracks{\sum_{\ell = 1}^{\infty}{1 \over \ell^{2}}
-\sum_{\ell = 1}^{\infty}{1 \over \pars{2\ell}^{2}}}
-{1 \over 4}\sum_{\ell = 1}^{\infty}{1 \over \ell^{2}}
=\half\ \underbrace{\sum_{\ell = 1}^{\infty}{1 \over \ell^{2}}}
_{\ds{\zeta\pars{2} = {\pi^{2} \over 6}}} =\ \color{#00f}{\large{\pi^{2} \over 12}}
\end{align}

$\zeta\pars{z}$ is the
Riemann Zeta Function .
A: Consider the contour integral
$$\oint_C dz \frac{z^2}{e^z+1}$$
where $C$ is the rectangle with vertices at $0$, $R$, $R+i 2 \pi$, and $i 2 \pi$, in that order, with a small semicircular indentation of radius $\epsilon$ at $z=i \pi$ into the rectangle.  Thus, the contour integral is equal to
$$\int_0^R dx \frac{x^2-(x+i 2 \pi)^2}{e^x+1} + i \int_0^{2 \pi} dy \frac{(R+i y)^2}{e^{R+i y}+1} \\+i PV \int_0^{2 \pi} dy \frac{y^2}{e^{i y}+1} + i \epsilon \int_{\pi/2}^{-\pi/2} d\phi \, e^{i \phi} \frac{(i \pi+\epsilon e^{i \phi})^2}{-e^{\epsilon e^{i \phi}}+1}$$ 
where $PV$ denotes the Cauchy principal value.  As $R \to \infty$, the second integral vanishes. As $\epsilon \to 0$, the fourth integral approaches $-i \pi^3$.  By Cauchy's theorem, the contour integral is zero.  Thus, we have
$$-i 4 \pi \int_0^{\infty} dx \frac{x}{e^x+1} + 4 \pi^2 \int_0^{\infty} \frac{dx}{e^x+1}\\+i PV \int_0^{2 \pi} dy \frac{y^2}{e^{i y}+1} - i \pi^3 = 0$$
Now,
$$i PV \int_0^{2 \pi} dy \frac{y^2}{e^{i y}+1} = \frac12 PV \int_0^{2 \pi} dy \, y^2 \,  \tan{\frac{y}{2}} + i \frac12 \int_0^{2 \pi} dy \, y^2 $$
Equating imaginary parts, we get that
$$-4 \pi \int_0^{\infty} dx \frac{x}{e^x+1} = \pi^3 - \frac{4 \pi^3}{3}$$
or
$$\int_0^{\infty} dx \frac{x}{e^x+1} = \frac{\pi^2}{12}$$
A: By way of diversity and since you admit complex variable methods note that what we have here is a special value of a Mellin transform:
$$\mathfrak{M}\left(\frac{1}{e^x+1}; s\right)
= \int_0^\infty \frac{1}{e^x+1} x^{s-1} dx.$$
To compute this transform we may proceed in an admittedly somewhat unorthodox manner by expanding the function being transformed into a series:
$$\int_0^\infty \frac{1}{e^x+1} x^{s-1} dx
= \int_0^\infty \frac{e^{-x}}{1+e^{-x}} x^{s-1} dx
= \int_0^\infty \sum_{q\ge 1} (-1)^{q+1} e^{-qx} x^{s-1} dx 
\\ = \sum_{q\ge 1} (-1)^{q+1}  \int_0^\infty  e^{-qx}  x^{s-1} dx
= \Gamma(s) \sum_{q\ge 1} \frac{(-1)^{q+1}}{q^s} 
= \Gamma(s) \left(1-\frac{2}{2^s}\right) \zeta(s).$$
Put $s=2$ to obtain
$$\Gamma(2) \times \frac{1}{2} \times \zeta(2) = \frac{\pi^2}{12}.$$
Here we have used the integral representation of the gamma function which I suppose is admissible in a complex variable method type answer.
A: We start from the contour integral evaluated through Cauchy's theorem: 
$$
\oint_{\Gamma_{\epsilon,R}}\mathrm{d}z\frac{z^2}{e^z+1} = 0
$$
where $\Gamma_{\epsilon, R}$ is the rectangle of vertices $0,\ R,\ R+i2\pi,\ i2\pi$, indented at the singularity $i\pi$ with a semicircle of radius $\epsilon$. 
Splitting it into its natural branches yields:
$$
\int_0^R\mathrm{d}x\frac{x^2}{e^x+1}+
\int_0^{2\pi}i\mathrm{d}y
\frac{(iy+R)^2}{e^{iy+R}+1}+
\int_R^0\mathrm{d}x\frac{(x+i2\pi)^2}{e^{x+i2\pi}+1}+
\left(\int_{2\pi}^{\pi+\epsilon}+
\int_{\pi-\epsilon}^0\right)i\mathrm{d}y\frac{(iy)^2}{e^{iy}+1}+
\int_{\pi/2}^{-\pi/2}i\epsilon e^{i\theta}\mathrm{d}\theta\frac{(\epsilon e^{i\theta}+i\pi)^2}{e^{i\pi + \epsilon e^{i\theta}}+1}=0.$$
Now we observe that the  first integral cancels out with the first term of the expansion of the third integral, furthermore the second one $\to0$ as $R\to \infty$. Our desired integral appears as the rectangular term of said expansion.
We are now left with the following calculations:
$$
-\int_{-\pi/2}^{\pi/2}i\epsilon e^{i\theta}\mathrm{d}\theta\frac{(\epsilon e^{i\theta}+i\pi)^2}{e^{i\pi + \epsilon e^{i\theta}}+1}=i\pi^2\int_{-\pi/2}^{\pi/2}\mathrm{d}\theta(-1+O(\epsilon)) \to_{\epsilon\to0}-i\pi^3
$$
and, using Euler's formula and splitting into real and imaginary part:
$$
i\int_0^{2\pi}\mathrm{d}y\frac{y^2}{\cos y+1+i\sin y}=i\int_0^{2\pi}\mathrm{d}y\ \frac{y^2}{2}- \int_0^{2\pi}\mathrm{d}y\frac{y^2\sin y}{2(\cos y +1}.
$$
Therefore, using the easy-to-verify $\int\mathrm{d}x\frac{1}{e^x+1}=\ln 2$:
$$
0-i4\pi I + 4\pi^2\ln2+i\int_0^{2\pi}\mathrm{d}y\ \frac{y^2}{2}-\int_0^{2\pi}\mathrm{d}y\frac{y^2\sin y}{2(\cos y +1}-i\pi^3=0.
$$
By comparison of real and imaginary part:
$$
I=\frac{\pi^2}{12}
$$
$$
\int_0^{2\pi}\mathrm{d}y\frac{y^2\sin y}{2(\cos y +1}=4\pi^2\ln 2
$$
