Integrating $\frac{x^3}{\exp(x)-1}$ from $0$ to $\infty$ While doing Physics and trying to prove the law of Stefan-Boltzmann from Plancks-law one comes to the integral 
\[ \int_0^\infty \frac{x^3}{\exp(x)-1}  \mathrm{d}x=\frac{\pi^4}{15} \]
and I would like to know how one calculates this value. One could try integration by parts to get rid of the $x^3$ part, but as Mathematica says the anti derivative is not very simple (there are some polylogarithms involved). Does anyone have an idea how to solve the Integral?
 A: Here, we develop an alternative, straightforward  approach that makes use of the Polylogarithm Functions.  To that end, we start with the integral 
$$I\equiv\int_0^{\infty}\frac{x^3}{e^x-1}dx$$
and enforce the substitution $x\to -\log (1-x)$ with $dx\to \frac{1}{1-x}dx$ to obtain
$$I=-\int_0^{1}\frac{(\log (1-x))^3}{x}dx$$
Next , we enter into a series of three integration by parts.
First, integration by parts with $u=(\log (1-x))^3$ and $v=\log x$ yields
$$I=-3\int_0^1 (\log (1-x))^2 \frac{\log x}{1-x}\,dx$$
Upon the second integration by parts with $u=(\log (1-x))^2$ and $v=\text{Li}_2(1-x)$, where $\text{Li}_2(x)\equiv\int_0^x \frac{\log (1-t)}{t}dt$ is the Dilogarithm Function, we obtain
$$I=-6\int_0^1 \log (1-x) \frac{\text{Li}_2(1-x)}{1-x}\,dx$$
And finally, a third and final integration by parts with $u=\log (1-x)$ and $v=\text{Li}_3(1-x)$, where $\text{Li}_3(x)\equiv\int_0^x \frac{\text{Li}_2 (1-t)}{t}dt$ is a Polylogarithm Function, reveals
$$\begin{align}
I&=6\int_0^1 \frac{\text{Li}_3(1-x)}{1-x}\,dx\\\\
&=-6\left.\text{Li}_4(1-x)\right|_0^1\\\\
&=6\text{Li}_4(1)\\\\
&=6\zeta(4)\\\\
&=\frac{\pi^4}{15}
\end{align}$$
as expected, where we have used The Relationship $\text{Li}_s(1)=\zeta(s)$, between the polylogarithm function and the Riemann-Zeta Function.
A: I will continue along Lee's comment to produce a full answer. 
$\displaystyle I=\int_0^\infty \frac{x^3}{e^{x}-1}dx$
$\displaystyle I=\int_0^\infty x^3(e^{-x}+e^{-2x}+\ldots)dx$
$\displaystyle I=\int_0^\infty x^3e^{-x}dx+\int_0^\infty x^3e^{-2x}dx+\ldots$
Now, 
$\displaystyle A=\int_0^\infty x^m e^{-nx} dx$
Let $x=\frac{u}{n},\ dx=\frac{1}{n}du$. Then, 
$\displaystyle A=\frac{1}{n^{m+1}}\int_0^\infty u^m e^{-u}du = \frac{m!}{n^{m+1}}  $
Thus, we get 
$\displaystyle I=3!\sum_{n=1}^\infty\frac{1}{n^{4}}$
$\displaystyle I=6\zeta(4)$
However, we know that $\displaystyle \zeta(4)=\frac{\pi^4}{90}$
Thus, $\displaystyle I=\frac{\pi^4}{15}$
