How to do integral $\int_{0}^{\infty} \frac{x^4 e^{x}}{(e^x-1)^2}\,dx$

The Riemann/Gamma way:$$\int_0^\infty\frac{x^4e^{-x}dx}{(1-e^{-x})^2}=\sum_{n\ge 1}n\int_0^\infty x^4e^{-nx}dx=\sum_{n\ge 1}\frac{4!}{n^4}=24\zeta(4)=\frac{4\pi^4}{15}.$$

Recall the integral representation of the Polylogarithm given by

$$\operatorname{Li}_s(z)~=~\frac1{\Gamma(s)}\int_0^\infty \frac{t^{s-1}}{e^t/z-1}\mathrm dt$$

Moreover note that

$$\frac{\partial}{\partial z}\left[\frac{t^{s-1}}{e^t/z-1}\right]=\frac{t^{s-1}e^t}{(e^t-z)^2}$$

which is precisely the structure of your integrand. Thus, by applying a variation of Feynman's Trick we obtain

\begin{align*} \frac{\mathrm d}{\mathrm dz}\Gamma(s)\operatorname{Li}_s(z)&=\frac{\mathrm d}{\mathrm dz}\int_0^\infty \frac{t^{s-1}}{e^t/z-1}\mathrm dt\\ \Gamma(s)\frac{\mathrm d}{\mathrm dz}\operatorname{Li}_s(z)&=\int_0^\infty \frac{\partial}{\partial z}\frac{t^{s-1}}{e^t/z-1}\mathrm dt\\ \Gamma(s)\frac{\operatorname{Li}_{s-1}(z)}z&=\int_0^\infty\frac{t^{s-1}e^t}{(e^t-z)^2}\mathrm dt \end{align*}

Now plugging in $s=5$ and $z=1$ yields to

\begin{align*} \Gamma(5)\frac{\operatorname{Li}_{5-1}(1)}1&=\int_0^\infty\frac{t^{5-1}e^t}{(e^t-1)^2}\mathrm dt\\ 4!\operatorname{Li}_4(1)&=\int_0^\infty\frac{t^{4}e^t}{(e^t-1)^2}\mathrm dt \end{align*}

Note that the RHS is the integral we are interested in whereas the LHS can be evaluated directly in terms of the Riemann Zeta Function. Overall this boils down to utilizing the value of $\zeta(4)$, as J.G.'s answer did aswell, since $\operatorname{Li}_s(1)=\zeta(s)$ as one can verify by using the series representation of the Polylogarithm.

$$\therefore~\int_0^\infty\frac{t^{4}e^t}{(e^t-1)^2}\mathrm dt~=~4!\operatorname{Li}_4(1)~=~\frac{4\pi^4}{15}$$