Find the value of $\int_{0}^{\infty}\frac{x^3}{e^x-1}\ln(e^x - 1)\,dx$

Mathematica says that the answer is $$\pi^2\zeta(3)+12\zeta(5)$$ I will try to figure out how this can be proven.

---

Added: Let me compute the 2nd integral in Ron Gordon's answer: \begin{align}\int_{0}^{\infty}\frac{x^3 e^{-x}}{1-e^{-x}}\ln(1-e^{-x})\,dx &=-\frac32\int_0^{\infty}x^2\ln^2(1-e^{-x})\,dx=\\&=-\frac32\left[\frac{\partial^2}{\partial s^2}\int_0^{\infty}e^{-sx}\ln^2(1-e^{-x})\,dx\right]_{s=0}=\\ &=-\frac32\left[\frac{\partial^2}{\partial s^2}\int_0^{1}t^{s-1}\ln^2(1-t)\,dt\right]_{s=0}=\\ &=-\frac32\left[\frac{\partial^4}{\partial s^2\partial u^2}\int_0^{1}t^{s-1}(1-t)^u\,dt\right]_{s=0,u=0}=\\ &=-\frac32\left[\frac{\partial^4}{\partial s^2\partial u^2}\frac{\Gamma(s)\Gamma(1+u)}{\Gamma(1+s+u)}\right]_{s=0,u=0}=\\ &=-\frac{1}{2}\left(\pi^2\psi^{(2)}(1)-\psi^{(4)}(1)\right). \end{align} To obtain the last expression, one should expand the ratio of gamma functions to 2nd order in $u$, then to expand the corresponding coefficient to 2nd order in $s$.

Then we can use that $\psi^{(2)}(1)=-2\zeta(3)$ and $\psi^{(4)}(1)=-24\zeta(5)$ (cf formula (15) here) to obtain the quoted result.

How about pulling factors of $e^{-x}$ from both the denominator and log terms? Then you end up with two separate integrals:

$$\int_0^{\infty}dx \frac{x^4 \, e^{-x}}{1-e^{-x}} + \int_0^{\infty}dx \frac{x^3 \, e^{-x}}{1-e^{-x}} \log{(1-e^{-x})}$$

In both cases, you Taylor expand the denominator in $e^{-x}$. For the first integral, this results in

$$\sum_{k=0}^{\infty} \int_0^{\infty}dx\, x^4 \, e^{-(k+1) x} = 4! \sum_{k=0}^{\infty} \frac{1}{(k+1)^5} = 24 \, \zeta(5) $$

For the second integral, you also need to Taylor expand the log term. This results in a double sum:

$$\begin{align}\sum_{k=0}^{\infty} \int_0^{\infty}dx\, x^3 \, e^{-(k+1) x} \log{(1-e^{-x})} &= -\sum_{k=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{m} \int_0^{\infty} dx \, x^3 e^{-(k+m) x}\\ &= - 3! \sum_{m=1}^{\infty} \frac{1}{m} \sum_{k=1}^{\infty} \frac{1}{(k+m)^4}\\ &= -\sum_{m=1}^{\infty} \frac{\psi^{(3)}(m+1)}{m} \end{align}$$

where $\psi$ is a polygamma function.

Using the change of variables $ u=e^{-x} $, we have

$$\int_{0}^{\infty}\frac{x^3}{e^x-1}\ln(e^x - 1)\,dx = \int _{0}^{1}\!{\frac { \left( \ln \left( u \right) \right) ^{3} \ln \left( 1-u \right) }{u-1 }}{du}- \int _{0}^{1}\!{\frac { \left( \ln \left( u \right) \right)^{4} }{u-1 }}{du}. $$

Now, just apply the technique which has been used to find the exact solution in this problem and the result will follow.