Conjecture ${\large\int}_0^\infty\left[\frac1{x^4}-\frac1{2x^3}+\frac1{12\,x^2}-\frac1{\left(e^x-1\right)x^3}\right]dx=\frac{\zeta(3)}{8\pi^2}$ I encountered the following integral and numerical approximations tentatively suggest that it might have a simple closed form:

$${\large\int}_0^\infty\left[\frac1{x^4}-\frac1{2x^3}+\frac1{12\,x^2}-\frac1{\left(e^x-1\right)x^3}\right]dx\stackrel{\color{gray}?}=\frac{\zeta(3)}{8\pi^2}\tag{$\diamond$}$$
  (Update: I fixed a typo: replaced $4\pi^2$ with $8\pi^2$ in the denominator)  

I have only about $800$ decimal digits that agree with the conjectured value, calculated using Mathematica. Unfortunately, its numerical algorithms become unstable when I try to increase precision. Maple refuses to numerically evaluate this integral altogether.
Obviously, the first three terms of the integrand have elementary antiderivatives, but I was not able to find a closed-form antiderivative (either elementary or using known special functions) for the last one.
I'm asking for your help in proving (or disproving) the $(\diamond)$.
 A: What about the Laplace transform? By using it, we have that our integral equals:
$$ I= \frac{1}{36}\int_{0}^{+\infty}\left(1-3 s+6 s^2-6 s^3 \psi'(1+s)\right)\,ds$$
and in this form Mathematica is perfectly able to state that $I=\frac{\zeta(3)}{\color{red}{8}\pi^2}$.
I just used:
$$\mathcal{L}^{-1}\left(\frac{1}{x^4}\right)=\frac{s^3}{6},\qquad \mathcal{L}\left(1-\frac{x}{2}+\frac{x^2}{12}-\frac{x}{e^x-1}\right)=\frac{1-3 s+6 s^2}{6 s^3}-\psi'(1+s) $$
together with:
$$ \int_{0}^{+\infty}f(x)g(x)\,dx = \int_{0}^{+\infty}(\mathcal{L} f)(s)(\mathcal{L}^{-1}g)(s)\,ds.$$
A: $$\text{for} \ \ \Re(s) > 1 : \qquad \int_0^\infty \frac{x^{s-1}}{e^x-1}dx = \Gamma(s) \zeta(s)$$
$\Gamma(s) \zeta(s)$ is meromorphic so that we can easily jump other its poles which are at $1,0,-2n+1$ for $n \in \mathbb{N}^*$.
the pole at $s=1$ is of residue $1$, so that :
$$\text{for} \ \ \Re(s) \in ]0;1[ : \qquad \int_0^\infty \frac{x^{s-1}}{e^x-1} - x^{s-2} dx = \Gamma(s) \zeta(s)$$
the pole at $s=0$ is of residue $\zeta(0) = -1/2$, so that :
$$\text{for} \ \ \Re(s) \in ]-1;0[ : \qquad \int_0^\infty \frac{x^{s-1}}{e^x-1} - x^{s-2} +\frac{x^{s-1}}{2} dx = \Gamma(s) \zeta(s)$$
the pole at $s=-1$ is of residue $-\zeta(-1) =  1/12$, so that :
$$\text{for} \ \ \Re(s) \in ]-3;-1[ : \qquad \int_0^\infty \frac{x^{s-1}}{e^x-1} - x^{s-2} + \frac{x^{s-1}}{2} - \frac{x^{s}}{12} dx = \Gamma(s) \zeta(s)$$
and finally when $s \to 0$ : $\Gamma(s-2) \approx \frac{1}{2s}$ and $\zeta(s-2) \approx s \zeta^{\prime}(-2) = -s\frac {2} {2 (2\pi)^{2}} \zeta (3)$ so that your integral is 
$$\lim_{s\to 0} - \Gamma(s-2) \zeta(s-2) = \frac { \zeta (3)} {8 \pi ^2}$$
https://en.wikipedia.org/wiki/Gamma_function#The_gamma_function_in_the_complex_plane
https://en.wikipedia.org/wiki/Particular_values_of_Riemann_zeta_function
https://fr.wikipedia.org/wiki/Fonction_z%C3%AAta_de_Riemann#Expression_int.C3.A9grale
A: From the summation identity of zeta function:
$$ \boxed{ {\,}\\ \quad \color{Blue}{\sum_{n=0}^{\infty}\frac{\Gamma(n+s)\zeta(n+s)}{(n+1)!}=0 \qquad\colon\space Re\{s\}\lt1} \quad \\{\,} } $$ 
$$ \begin{align} 
\Gamma(s)\zeta(s) &= -\sum_{n=1}^{\infty}\frac{\Gamma(n+s)\zeta(n+s)}{(n+1)!} = -\int_{0}^{\infty}\frac{x^{s-2}}{e^x-1}\left(\sum_{n=1}^{\infty}\frac{x^{n+1}}{(n+1)!}\right)\,dx \\[2mm] 
&= -\int_{0}^{\infty}\frac{x^{s-2}}{e^x-1}\left(e^x-1+x\right)\,dx = \int_{0}^{\infty}x^{s-2}\left(\frac{x}{e^x-1}-1\right)\,dx \\[4mm] 
\Gamma(s-1)\zeta(s-1) &= -\frac{\Gamma(s)\zeta(s)}{2!}-\sum_{n=2}^{\infty}\frac{\Gamma(n+s)\zeta(n+s)}{(n+1)!} \\[2mm] 
&= \int_{0}^{\infty}x^{s-3}\left(\frac{x}{e^x-1}-1+\frac{x}{2}\right)\,dx \qquad\cdots\,\implies 
\end{align} $$ 

$$ \color{blue}{\Gamma(s-N)\zeta(s-N)=\int_{0}^{\infty}x^{s-N-2}\left[\frac{x}{e^x-1}-\left(\sum_{n=0}^{N}B_{n}\frac{x^n}{n!}\right)\right]\,dx} $$ 
$$ {\small \,0\lt\,Re\{s\}\,\lt1 ,\quad N\in\{\,0,\,1,\,2,\,\cdots\,\} ,\quad B_{n}\,\,{Bernoulli\,Number} ,\quad B_{1}=-1/2} $$ 

 
$$ \begin{align} 
\color{red}{I} &= \int_{0}^{\infty}\left[\frac{1}{x^4}-\frac{1}{2\,x^3}+\frac{1}{12\,x^2}-\frac{1}{\left(e^x-1\right)\,x^3}\right]\,dx \\[3mm] 
&= -\int_{0}^{\infty}x^{-4}\left[\frac{x}{e^x-1}-1+\frac{x}{2}-\frac{x^2}{12}\right]\,dx \\[3mm] 
&= -\int_{0}^{\infty}x^{\color{red}{0-2}-2}\left[\frac{x}{e^x-1}-\left(1\frac{x^0}{0!}-\frac{1}{2}\frac{x^1}{1!}+\frac{1}{6}\frac{x^2}{2!}\right)\right]\,dx \\[3mm] 
&= -\lim_{s\rightarrow0}\Gamma(s-2)\zeta(s-2)=-\frac{\zeta'(-2)}{2}=\color{red}{\frac{\zeta(3)}{8\pi^2}} 
\end{align} $$ 
