Integration evaluate $\int_0^\infty\frac{x^2}{e^x-1}\mathrm{dx}$. can I evaluate this integration ?

$$\int_0^\infty\dfrac{x^2}{e^x-1}\mathrm{dx}$$

how?
thank you
 A: \begin{align}
\int_0^\infty\dfrac{x^2}{e^x-1}\mathrm{dx} 
&= \int_0^\infty\dfrac{x^2e^{-x}}{1-e^{-x}}\mathrm{dx} \\
&= \int_0^\infty\mathrm{dx}~ x^2e^{-x}\sum_{n=0}^{\infty}e^{-nx}  \\
&= \sum_{n=0}^{\infty}\int_0^\infty x^2e^{-(n+1)x}\mathrm{dx}  \\
&= \sum_{n=0}^{\infty}\dfrac{1}{(n+1)^3}\Gamma(3) \\
&= \color{blue}{2\zeta(3)}
\end{align}
A: Here is an approach which avoids series and is based on the polylogarithm function.
We begin by rewriting the integral as
$$I = \int_0^\infty \frac{x^2}{e^x - 1} \, dx = \int^\infty_0 \frac{x^2 e^{-x}}{1 - e^{-x}} \, dx.$$
Now as the polylogarithm function of order zero is given by
$$\text{Li}_0 (x) = \frac{x}{1 - x},$$
if $x$ is replaced with $e^{-x}$ in the above expression, we have
$$\text{Li}_0 (e^{-x}) = \frac{e^{-x}}{1 - e^{-x}},$$
and the integral can be rewritten in terms of the polylogarithm function as
$$I = \int^\infty_0 x^2 \text{Li}_0 (e^{-x}) \, dx. \tag1$$
As the derivative for the polylogarithmic function of order $s$ is given by
$$\frac{d}{dx} \text{Li}_s (x) = \frac{\text{Li}_{s - 1} (x)}{x},$$
one has
$$\frac{d}{dx} \text{Li}_s (e^{-x}) = - \text{Li}_{s - 1} (e^{-x}),$$
and we immediately see that
$$\int \text{Li}_s (e^{-x}) \, dx = - \text{Li}_{s + 1} (e^{-x}) + C. \tag2$$
We now make use of (2) and repeatedly integrate (1) by parts. Doing so yields
\begin{align*}
I &= 2 \int^\infty_0 x \text{Li}_1 (e^{-x}) \, dx \quad \text{(by parts)}\\
&=  2 \int^\infty_0 \text{Li}_2 (e^{-x}) \, dx \quad \text{(by parts again)}\\
&= 2 \, \text{Li}_3 (1).
\end{align*}
Making use of the following result
$$\text{Li}_s (1) = \zeta (s), \quad s > 1,$$
where $\zeta(x)$ is the Riemann zeta function we finally obtain
$$\int^\infty_0 \frac{x^2}{e^x - 1} \, dx = 2 \zeta (3).$$
