= Solution
For $x>0$, differentiate the absolutely convergent <geometric series> to obtain
$$
(1-e^{-x})^{-2}=\sum_{j=0}^\infty(j+1)e^{-jx},\qquad
(e^x-1)^{-2}=\sum_{m=2}^\infty(m-1)e^{-mx}.
$$
If $\sigma=\Re z>1$, the sum of <integrals> of absolute values is
$$
\sum_{m\ge2}(m-1)\int_0^\infty x^\sigma e^{-mx}\,dx
=\Gamma(\sigma+1)\sum_{m\ge2}\frac{m-1}{m^{\sigma+1}}<\infty.
$$
Thus <Fubini theorem> permits termwise integration, even for complex $z$. The substitution $u=mx$ and the definition of the <Gamma function> give
$$
\boxed{\int_0^\infty\frac{x^z}{(e^x-1)^2}\,dx
=\Gamma(z+1)\sum_{m\ge2}\left(m^{-z}-m^{-z-1}\right)
=\Gamma(z+1)\bigl(\zeta(z)-\zeta(z+1)\bigr).}
$$
The missing $m=1$ terms cancel in the difference of the two <Riemann zeta functions>.
Back to article page
