I'm trying to figure out how to evaluate the following: $$ J=\int_{0}^{\infty}\frac{x^3}{e^x-1}\ln(e^x - 1)\,dx $$ I'm tried considering $I(s) = \int_{0}^{\infty}\frac{x^3}{(e^x-1)^s}\,dx\implies J=-I'(1)$, but I couldn't figure out what $I(s)$ was. My other idea was contour integration, but I'm not sure how to deal with the logarithm. Mathematica says that $J\approx24.307$.
I've asked a similar question and the answer involved $\zeta(s)$ so I suspect that this one will as well.
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.