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?
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?
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}$$
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)$$
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] $$
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}