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)$.

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.$$

$$\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