Prove $\int^\infty_0\frac x{e^x-1}dx=\frac{\pi^2}{6}$

The standard way of attacking this integral is to Taylor expand the denominator:

$$\int_0^{\infty} dx \frac{x^{u-1}}{e^x-1} = \sum_{k=0}^{\infty} \int_0^{\infty} dx\, x^{u-1} \, e^{-(k+1) x} = \Gamma(u) \sum_{k=0}^{\infty} \frac{1}{(k+1)^u} = \Gamma(u) \zeta(u)$$

I really don't know what insight you would gleam by using a contour integral in the complex plane. You could use a keyhole contour, and the you would get residues at integers along the imaginary axis so that you will recover the result above. But the method above is far easier and more intuitive.

UPDATED
You may indeed use contour integration as nicely exposed in Ablowitz and Fokas "Complex Variables" p. $244$.

The idea is to evaluate $\;\displaystyle I(p):=PV \int_{-\infty}^\infty \frac {e^{px}}{1-e^x}dx,\quad(0<p<1)\quad$ using the contour :

The integrals at the right and left side vanish as $R\to +\infty$ and the four bottom and top horizontal lines will add the contribution : $$\left(\int_{-R}^{-\epsilon_1}+\int_{\epsilon_1}^R\right)\frac{e^{px}}{1-e^x}dx+e^{2\pi ip}\left(\int_{R}^{-\epsilon_2}+\int_{\epsilon_2}^{-R}\right)\frac{e^{px}}{1-e^x}dx$$ that is (as $\;R\to +\infty$ and $\epsilon_i\to 0$) : $$\tag{1}\left(1-e^{2\pi ip}\right)PV \int_{-\infty}^\infty \frac{e^{px}}{1-e^x}dx$$while the semi-circular contours will give, as $\epsilon_1 \to 0^+$ and $\epsilon_2 \to 0^+$ :

\begin{align} \tag{2}\lim_{\epsilon_1\to 0}\int_{C\epsilon_1}\frac{e^{pz}}{1-e^z}dz&=-i\pi\left(\frac{e^{pz}}{-e^z}\right)_{z=0}=i\pi\\ \tag{3}\lim_{\epsilon_2\to 0}\int_{C\epsilon_2}\frac{e^{pz}}{1-e^z}dz&=-i\pi\left(\frac{e^{pz}}{-e^z}\right)_{z=2\pi i}=i\pi\,e^{2\pi i z}\\ \end{align}

Since no singularity appears in its inside the whole contour will have an integral of value $0$ equal to the sum of $(1),\;(2)$ and $(3)$. This gives : $$I(p)=PV \int_{-\infty}^\infty \frac {e^{px}}{1-e^x}dx=-i\pi \frac{1+\,e^{2\pi i z}}{1-\,e^{2\pi i z}}$$ or simply : $$\tag{4} \boxed{\displaystyle I(p)=\pi \cot( \pi p)}$$ Now $(4)$ multiplied by $p$ should give Euler's generating function for $\zeta($even$)$ : $$\tag{5}\displaystyle\pi\;p\;\cot(\pi\;p)=1-2\sum_{n=1}^\infty \zeta(2n)\;p^{2n}$$ which allows to find $\zeta(2n)$ just by using the expansion of $\pi\,x\;\cot(\pi\,x)$ : $$\tag{6}\pi\,x\;\cot(\pi\,x)=1-2\left(\frac {\pi^2}6x^2+\frac{\pi^4}{90}x^4+\frac{\pi^6}{945}x^6+\cdots\right)$$

Since we had an integral from $0$ to $+\infty$ and not a Cauchy principal value from $-\infty$ to $+\infty$ it will be convenient to rewrite $I(p)$ as :

\begin{align} \tag{7}I(p)-\frac 1p&=\pi \cot( \pi p)-\frac 1p=-2\sum_{n=1}^\infty \zeta(2n)\;p^{2n-1}\\ &=PV \int_{-\infty}^\infty \frac {e^{px}}{1-e^x}dx-\int_0^\infty e^{-px}dx\\ &=\lim_{\epsilon\to 0} \int_{\epsilon}^\infty \frac {e^{px}}{1-e^x}-e^{-px}dx+\int_{-\infty}^{-\epsilon} \frac {e^{py}}{1-e^y}dy\\ (\text{setting}\;\;y:=-x)\;\\ &=\lim_{\epsilon\to 0} \int_{\epsilon}^\infty \frac {e^{px}}{1-e^x}-e^{-px}+\frac {e^{-px}e^x}{(1-e^{-x})e^x}dx\\ &=\lim_{\epsilon\to 0} \int_{\epsilon}^\infty \frac {e^{px}}{1-e^x}+\frac {e^{-px}e^x-e^{-px}(e^x-1)}{e^x-1}dx\\ \end{align} with the neat result :
$$\tag{8}\boxed{\displaystyle I(p)-\frac 1p= \int_0^\infty \frac {e^{px}-e^{-px}}{1-e^x}dx}$$

This integral may be used to obtain the coefficients of $p^{2n-1}$ in $(7)$ since : $$\left(\frac d{dp}\right)^{2n-1} \left(I(p)-\frac 1p\right)=\int_0^\infty \frac {x^{2n-1}\left(e^{px}+e^{-px}\right)}{1-e^x}dx$$ with the limit as $p\to 0$ given by : $$\int_0^\infty \frac {x^{2n-1}\;2}{1-e^x}dx=-2\,(2n-1)!\,\zeta(2n)$$ corresponding to the classical formula for $n$ positive integer : $$\tag{9}\zeta(2n)\Gamma(2n)=\int^\infty_0\frac{x^{2n-1}}{e^x-1}dx$$

$\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}}$ \begin{align} &\color{#00f}{\large\int_{0}^{\infty}{x \over \expo{x} - 1}\,\dd x}= \int_{x\ =\ 0}^{x\ \to\ \infty}x\,\dd\bracks{\ln\pars{1 - \expo{-x}}} =-\int_{0}^{\infty}\ln\pars{1 - \expo{-x}}\,\dd x \\[3mm]&=\sum_{n = 1}^{\infty}{1 \over n}\int_{0}^{\infty}\expo{-nx}\,\dd x =\sum_{n = 1}^{\infty}{1 \over n^{2}}\ \overbrace{\int_{0}^{\infty}\expo{-x}\,\dd x}^{\ds{=\ 1}}\ =\sum_{n = 1}^{\infty}{1 \over n^{2}}=\color{#00f}{\large{\pi^{2} \over 6}} \end{align}

Let $$f(z) = \frac{z}{e^z-1}$$ $f$ has poles at $z = 2 \pi i n$ where $n \in \mathbb Z$ and $n \neq 0$. Let $C$ be the rectangular contour with vertices at $0, R, R+i\pi$ and $i\pi$ but indented at 0. Because the residue is $0$ there, the indentation does not add to the value of the integral.

Then, as $R \to \infty$:

$$0=\oint_C f(z)\, dz = \int_{0}^\infty f(z)\, dz-\int_0^\infty f(z+\pi i)\,dz-i\int_0^\pi f(iz)\,dz+\pi^2=\\ \int_{0}^\infty f(z)\, dz+\int_0^\infty \frac{z+\pi i}{e^z+1}\,dz+\int_0^\pi \frac{z}{e^{iz}-1}\,dz $$

Taking the real parts of both sides yeilds

$$0 = \int_{0}^\infty f(z)\, dz+\int_{0}^\infty \frac{z}{e^z+1}\, dz-\int_0^\pi \frac{x}{2}\,dz$$

so

$$\frac{\pi^2}{4}=\int_{0}^\infty f(z)\, dz+\int_{0}^\infty \frac{z}{e^z+1}\, dz$$

By summing and integrating both series as was done by Ron Gordon, we find

$$I=\int_{0}^\infty f(z)\, dz=\sum_{n=1}^\infty \frac{1}{n^2}$$ $$J=\int_{0}^\infty f(z)\, dz=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n^2} = \frac{I}{2}$$

and therefore

$$\frac{\pi^2}{4}=I+\frac{I}{2}$$

therefore

$$I=\int_{0}^\infty f(z)\, dz = \frac{\pi^2}{6}$$