Evaluate $\int_0^{\infty}\frac{x^2}{(e^x-1)^2}dx$

I want to do

$$\int_0^\infty\!\!\! \frac{x^2}{(e^x-1)^2}\text{d}x$$

I know how to do $\int_0^\infty \frac{x}{(e^x-1)}\text{d}x$ by using $-\frac{\ln(1-x)}{x} = \sum_{n=1}^\infty \frac{x^{n-1}}{n}$. But when I try to do the same trick I get a double sum which I have no idea of how to simply.

The result should be able to be presented as an expression involving $\zeta(3)$.

One more approach: the integral is $$\int_0^\infty\frac{x^2 e^{-2x}}{(1-e^{-x})^2}dx=\sum_{n=2}^\infty(n-1)\int_0^\infty x^2 e^{-nx}dx=2\sum_{n=2}^\infty (n^{-2}-n^{-3})=2(\zeta(2)-\zeta(3)).$$

Lets start with the integral representation of the Polylogarithm $\operatorname{Li}_s(z)$ which is given by

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

furthermore note the relations

$$\frac{\mathrm d}{\mathrm dz}\operatorname{Li}_s(z)=\frac{\operatorname{Li}_{s-1}(z)}{z}\text{ and }\operatorname{Li}_s(1)=\zeta(s)\tag2$$

Differentiate $(1)$ w.r.t. $z$ by using the first property of $(2)$ yields

$$\begin{align} \frac{\mathrm d}{\mathrm dz}\Gamma(s)\operatorname{Li}_s(z)&=\frac {\mathrm d}{\mathrm dz}\int_0^{\infty}t^{s-1}\frac{z}{e^t-z}\mathrm dt\\ \Gamma(s)\frac{\operatorname{Li}_{s-1}(z)}{z}&=\int_0^{\infty}t^{s-1}\frac{(e^t-z)+z}{(e^t-z)^2}\mathrm dt\\ \Gamma(s)\frac{\operatorname{Li}_{s-1}(z)}{z}&=\int_0^{\infty}\frac{t^{s-1}}{e^t-z}\mathrm dt+z\int_0^{\infty}\frac{t^{s-1}}{(e^t-z)^2}\mathrm dt\\ \Gamma(s)\frac{\operatorname{Li}_{s-1}(z)}{z}&=\Gamma(s)\frac{\operatorname{Li}_{s}(z)}{z}+z\int_0^{\infty}\frac{t^{s-1}}{(e^t-z)^2}\mathrm dt\\ \therefore~\int_0^{\infty}\frac{t^{s-1}}{(e^t-z)^2}\mathrm dt&=\Gamma(s)\left[\frac{\operatorname{Li}_{s-1}(z)}{z^2}-\frac{\operatorname{Li}_{s}(z)}{z^2}\right] \end{align}$$

Setting $s=3$ and $z=1$ and by including the second relation of $(2)$ gives us

$$\begin{align} \int_0^{\infty}\frac{t^{3-1}}{(e^t-1)^2}\mathrm dt&=\Gamma(3)\left[\frac{\operatorname{Li}_{3-1}(1)}{1^2}-\frac{\operatorname{Li}_{3}(1)}{1^2}\right]\\ \int_0^{\infty}\frac{t^{2}}{(e^t-1)^2}\mathrm dt&=2[\zeta(2)-\zeta(3)] \end{align}$$

$$\therefore~\int_0^{\infty}\frac{t^{2}}{(e^t-1)^2}dt=\frac{\pi^2}3-2\zeta(3)=0.885~754...$$

Where the integral on the left equals your integral by setting $t=x$. Numerically WolframAlpha gives the same solution.

Your approach to $\int_0^{\infty}\frac x{e^x-1}dx$ is also given through the special case $s=2$. Note that the derivative of the Dilogarithm $\operatorname{Li}_2(z)$ is given by the term $-\ln(1-x)/x$. Therefore your attempt was not wrong at all.

Withouth invoking the polylogarithm, for any $n\in\mathbb{N}$ and any $a>0$ we have

$$\begin{eqnarray*} I(n,a)=\int_{0}^{+\infty}\frac{x^n}{e^{ax}-1}\,dx&=&\sum_{m\geq 1}\int_{0}^{+\infty}x^n e^{-ma x}\,dx\\&=&\sum_{m\geq 1}\frac{n!}{(ma)^{n+1}}=\frac{n!}{a^{n+1}}\zeta(n+1)\end{eqnarray*}$$ and by differentiation under the integral sign $$ J(n,a)=\frac{\partial}{\partial a} I(n,a)=\int_{0}^{+\infty}\frac{-x^{n+1} e^{ax}}{(e^{ax}-1)^2}\,dx = -\frac{(n+1)!}{a^{n+2}}\zeta(n+1)$$ such that $$\begin{eqnarray*} \int_{0}^{+\infty}\frac{x^{n+1}}{(e^{ax}-1)^2}\,dx &=& -\left[I(n+1,a)+J(n,a)\right]\\&=&\frac{(n+1)!}{a^{n+1}}\zeta(n+1)-\frac{(n+1)!}{a^{n+2}}\zeta(n+2)\end{eqnarray*}$$ and by setting $a=n=1$ $$ \int_{0}^{+\infty}\frac{x^{2}}{(e^{x}-1)^2}\,dx =\color{red}{2\left(\zeta(2)-\zeta(3)\right)}=\sum_{n\geq 1}\frac{2n}{(n+1)^3}=\sum_{n\geq 1}\frac{6n+5(-1)^{n}}{n^3\binom{2n}{n}}. $$