Evaluate the integral $\int_{0}^{+\infty}\left(\frac{x^2}{e^x-1}\right)^2dx$ 
Find the value of the integral $$\int\limits_{0}^{+\infty}\left(\frac{x^2}{e^x-1}\right)^2\;\mathrm{d}x\;.$$

We can let $x=\ln{t}$ to get
$$\int\limits_{1}^{+\infty}\frac{(\ln{t})^4}{t(t-1)^2}\;\mathrm{d}t\;.$$
But how can we evaluate it from there?
 A: Since
$$ \frac{1}{(e^x-1)^2} = \sum_{n\geq 1} n e^{-(n+1)x} \tag{1}$$
and 
$$ \int_{0}^{+\infty} nx^4 e^{-(n+1)x}\,dx = \frac{24\,n}{(1+n)^5}\tag{2}$$
it follows that:

$$ \int_{0}^{+\infty}\left(\frac{x^2}{e^{x}-1}\right)^2\,dx = 24\sum_{n\geq 1}\frac{n}{(1+n)^5} = \color{red}{\frac{4}{15}\left(\pi^4-90\,\zeta(5)\right)}.\tag{3}$$

A: Using your substitution, we can continue in the following way:
$$I=\int_{1}^{+\infty}\dfrac{\ln^4{t}}{t(t-1)^2}dt$$
Let's make another substitution:
$$y=\frac{1}{t}$$
$$I=\int_{0}^{1}\frac{y \ln^4{y}}{(1-y)^2}dy$$
Integrating by parts with:
$$u=y \ln^4{y},~~~~~~~dv=\frac{dy}{(1-y)^2}$$
We obtain:
$$I=\int_{0}^{1}\frac{\ln^4{y}}{1-y}dy+4\int_{0}^{1}\frac{\ln^3{y}}{1-y}dy$$
I leave the proof of $u(0)v(0)=u(1)v(1)=0$ out for now.
Using the general formula:
$$\int^1_0 \frac{\log^n x}{1-x}dx=(-1)^n~ n!~ \zeta(n+1)$$
We get:

$$I=24 (\zeta(4)- \zeta(5))=\frac{4 \pi^4}{15}-24 \zeta(5)$$

A: The integral representation of polylogarithm function is
$$
\mathrm{Li}_\nu(z)=\frac{z}{\Gamma(\nu)}\int_0^\infty\frac{x^{\nu-1}}{\mathrm e^x-z}\mathrm d x
$$
and differentiating with respect to $z$ we have
$$
\mathrm{Li}'_\nu(z)=\frac{1}{\Gamma(\nu)}\int_0^\infty\frac{x^{\nu-1}}{\mathrm e^x-z}\mathrm d x+\frac{z}{\Gamma(\nu)}\int_0^\infty\frac{x^{\nu-1}}{(\mathrm e^x-z)^2}\mathrm d x
$$
thus
$$
\mathrm{Li}'_{\nu+1}(1)=\mathrm{Li}_{\nu+1}(1)+\frac{1}{\Gamma(\nu+1)}\int_0^\infty\frac{x^{\nu}}{(\mathrm e^x-1)^2}\mathrm d x
$$
and observing that $\mathrm{Li}'_{\nu+1}(z)=\frac{1}{z}\mathrm{Li}_{\nu}(z)$ and $\mathrm{Li}_{\nu}(1)=\zeta(\nu)$ we have
$$
\int_0^\infty\frac{x^{\nu}}{(\mathrm e^x-1)^2}\mathrm d x=\Gamma(\nu+1)\Big[\mathrm{Li}_{\nu}(1)-\mathrm{Li}_{\nu+1}(1)\Big]=\Gamma(\nu+1)\Big[\zeta(\nu)-\zeta(\nu+1)\Big]$$
and for $\nu=4$
$$
\int_0^\infty\frac{x^{4}}{(\mathrm e^x-1)^2}\mathrm d x=\Gamma(5)\Big[\zeta(4)-\zeta(5)\Big]=24\left[\frac{\pi^2}{90}-\zeta(5)\right]=\frac{4}{15}\left[\pi^2-90\,\zeta(5)\right]
$$
A: Let
\begin{equation}
f(x) = \frac{1}{(\mathrm{e}^{x}-1)^{2}}
\end{equation}
Then the Mellin transform of the function is
\begin{equation}
\mathcal{M}[f](s) = \int\limits_{0}^{\infty} \frac{x^{s-1}}{(\mathrm{e}^{x}-1)^{2}} \mathrm{d} x
= \Gamma(s)[\zeta(s-1) - \zeta(s)]
\end{equation}
For $s=5$, we have
\begin{equation}
\int\limits_{0}^{\infty} \frac{x^{4}}{(\mathrm{e}^{x}-1)^{2}} \mathrm{d} x = \Gamma(5)[\zeta(4) - \zeta(5)] 
= \frac{4\pi^{4}}{15} - 24\zeta(5)
\end{equation}
