Evaluation of $\int^{\infty}_{0}\frac{x^3}{e^x-1}dx$ 
Evaluation of $$\int^{\infty}_{0}\frac{x^3}{e^x-1}dx$$

$\bf{My\; Try::}$ Let $\displaystyle I = \int^{\infty}_{0}\frac{x^3}{e^x-1}dx$
Now put $\displaystyle e^x=\frac{1}{t}\;,$ Then $\displaystyle e^xdx = -\frac{1}{t^2}dt$ and $x=-\ln t$ and changing limits, We get
$$I = \int^{0}_{1}\frac{\ln^3 t}{1-t}dt = -\int^{1}_{0}\frac{\ln^3 t}{1-t}dt=-\int^{1}_{0}\ln^3(t)\sum^{\infty}_{n=0}t^ndt=\sum^{\infty}_{n=0}\int^{1}_{0}\ln^3 (t)t^ndt$$
Now Using By parts, We get $$I = 0+3\sum^{\infty}_{n=0}\frac{1}{n+1}\int^{1}_{0}\ln^2 (t)\cdot t^ndt$$
Again using By parts, We get $$I = 0-6\sum^{\infty}_{n=0}\frac{1}{(n+1)^2}\int^{1}_{0}\ln(t)\cdot t^ndt$$
Againg using by parts, We get $$I = 6\sum^{\infty}_{n=0}\frac{1}{(n+1)^3}\int^{1}_{0}t^ndt = \sum^{\infty}_{n=0}\frac{6}{(n+1)^4}=6\zeta(4)$$
Now how can i solve it, Help required, Thanks
 A: Making the problem more general, consider $$I=\int^{\infty}_{0}\frac{x^n}{e^{ax}-1}dx\qquad (a>0)$$ Changing variable $ax=t$ leads to $$I=\frac 1 {a^{n+1}}\int^{\infty}_{0}\frac{t^n}{e^{t}-1}dx$$
Using what  mathlove answered here, we then end with $$I=\frac{\zeta (n+1)\, \Gamma (n+1)}{a^{n+1}}$$ and nice expressions when $n$ is odd.
A: Expanding Claude Leibovici's solution,
\begin{align}
\int\limits_{0}^{\infty} \frac{x^{n}}{\mathrm{e}^{ax}-1} dx
&= \int\limits_{0}^{\infty} \frac{\mathrm{e}^{-ax} x^{n}}{1-\mathrm{e}^{-ax}} dx \\
&= \int\limits_{0}^{\infty} \mathrm{e}^{-ax} x^{n} \sum\limits_{k = 0}^{\infty} \mathrm{e}^{-kax} dx \\
&= \sum\limits_{k = 0}^{\infty} \int\limits_{0}^{\infty} x^{n} \mathrm{e}^{-(a+ak)x} dx
\end{align}
To evaluate the integral, let $(a+ak)x = z$
\begin{equation}
\int\limits_{0}^{\infty} x^{n} \mathrm{e}^{-(a+ak)x} dx
= \frac{1}{a^{n+1}} \frac{1}{(1+k)^{n+1}} \int\limits_{0}^{\infty} z^{n} \mathrm{e}^{-z} dz
= \frac{1}{a^{n+1}} \frac{1}{(1+k)^{n+1}} \Gamma(n+1)
\end{equation}
Now we have
\begin{equation}
\int\limits_{0}^{\infty} \frac{x^{n}}{\mathrm{e}^{ax}-1} dx
= \frac{1}{a^{n+1}} \Gamma(n+1) \sum\limits_{k = 1}^{\infty} \frac{1}{k^{n+1}}
= \frac{1}{a^{n+1}} \Gamma(n+1) \zeta(n+1)
\end{equation}
The original problem is
\begin{equation}
\int\limits_{0}^{\infty} \frac{x^{3}}{\mathrm{e}^{x}-1} dx
= \Gamma(4) \zeta(4) = 6\frac{\pi^{4}}{90} = \frac{\pi^{4}}{15}
\end{equation}
