Definite integral involving polylogarithm
I wonder if there is not a problem somewhere.
I agree that $$\sum_{k=1}^\infty \frac{H_k}{k^4}=-\zeta(2)\zeta(3)+\text{something}$$ Setting $$\text{something}=\int_0^1 \frac{\text{Li}_2^2(x)}{x}dx$$ may be not correct since the numerical integration leads to $0.843825$ while, numerically $$\sum_{k=1}^\infty \frac{H_k}{k^4}=1.13348$$ making $$\text{something}=3.11078$$ which is identified as $3\zeta(5)$.
It could be good that you describe the intermediate steps.
The explicit values of all the Euler sums $\sum_{n\geq 1}\frac{H_n}{n^s}$ are well known and related to convolutions of $\zeta$ values. Flajolet and Salvy showed how to derive them from the residue theorem, for instance.
Here I will outline a different technique, less efficient but more elementary. Let us assume that $a,b$ are distinct positive real numbers and recall that $\sum_{n\geq 1}H_n z^n = -\frac{\log(1-z)}{1-z}$. By partial fraction decomposition $$ \frac{1}{(n+a)(n+b)}=\frac{\frac{1}{b-a}}{n+a}+\frac{\frac{1}{a-b}}{n+b}$$ hence $$ \sum_{n\geq 1}\frac{H_n}{(n+a)(n+b)} = -\int_{0}^{1}\frac{\log(1-x)}{1-x}\left[\frac{x^{a-1}}{b-a}+\frac{x^{b-1}}{a-b}\right]\,dx$$ where the resulting integral can be computed by integration by parts and by differentiating Euler's Beta function. The outcome is: $$\sum_{n\geq 1}\frac{H_n}{(n+a)(n+b)} = \frac{H_{a-1}^2-H_{b-1}^2-\psi'(a)-\psi'(b)}{2(a-b)}$$ hence by considering the limit as $b\to a$: $$ \sum_{n\geq 1}\frac{H_n}{(n+a)^2}=H_{a-1}\psi'(a)-\tfrac{1}{2}\psi''(a). $$ In order to compute $\sum_{n\geq 1}\frac{H_n}{n^s}$, it is enough to apply $\left(\lim_{a\to 0^+}\frac{d^{s-2}}{da^{s-2}}\right)$ to both sides of the previous line. In the $s=4$ case we get: $$ \sum_{n\geq 1}\frac{H_n}{n^4}=\frac{1}{6}\lim_{a\to 0^+}\left[3\,\psi'(a)\,\psi''(a)+H_{a-1}\,\psi'''(a)-\tfrac{1}{2}\psi^{IV}(a)\right]=\color{blue}{3\,\zeta(5)-\zeta(2)\,\zeta(3)}. $$
$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} \sum_{k = 1}^{\infty}{H_{k} \over k^{4}} & = \sum_{k = 1}^{\infty}H_{k}\ \overbrace{\bracks{-\,{1 \over 6}\int_{0}^{1}\ln^{3}\pars{x}\,x^{k - 1}\,\dd x}} ^{\ds{1 \over k^{4}}} \\[5mm] & = -\,{1 \over 6}\int_{0}^{1}\ln^{3}\pars{x}\ \overbrace{\sum_{k = 1}^{\infty}H_{k}\,x^{k}} ^{\ds{-\,{\ln\pars{1 - x} \over 1 - x}}}\ \,{\dd x \over x} \\[5mm] & = {1 \over 6}\int_{0}^{1}{\ln^{3}\pars{x}\ln\pars{1 - x} \over x\pars{1 - x}} \,\dd x \\[5mm] & = {1 \over 6}\int_{0}^{1}{\ln^{3}\pars{x}\ln\pars{1 - x} \over x}\,\dd x + {1 \over 6}\int_{0}^{1}{\ln^{3}\pars{x}\ln\pars{1 - x} \over 1 - x}\,\dd x \\[5mm] & = {1 \over 6}\ \underbrace{\int_{0}^{1}{\ln^{3}\pars{x}\ln\pars{1 - x} \over x}\,\dd x} _{\ds{\mc{I}_{1}}}\ +\ {1 \over 6}\ \underbrace{\int_{0}^{1}{\ln^{3}\pars{1 - x}\ln\pars{x} \over x}\,\dd x} _{\ds{\mc{I}_{2}}} \\[5mm] & = {\mc{I}_{1} + \mc{I}_{2} \over 6}\label{1}\tag{1} \end{align}
$\ds{\Large \mc{I}_{1} = ?}$. \begin{align} \mc{I}_{1} & = -\int_{0}^{1}\mrm{Li}_{2}'\pars{x}\ln^{3}\pars{x}\,\dd x = 3\int_{0}^{1}\mrm{Li}_{3}'\pars{x}\ln^{2}\pars{x}\,\dd x = -6\int_{0}^{1}\mrm{Li}_{4}'\pars{x}\ln\pars{x}\,\dd x \\[5mm] & = 6\int_{0}^{1}\mrm{Li}_{5}'\pars{x}\,\dd x = 6\,\mrm{Li}_{5}\pars{1} \implies \bbx{\large \mc{I}_{1} = 6\,\zeta\pars{5}}\label{2}\tag{2} \end{align}
$\ds{\Large \mc{I}_{2} = ?}$. \begin{align} \mc{I}_{2} & = \left.\partiald[3]{}{\mu}\partiald{}{\nu} \int_{0}^{1}{\bracks{\pars{1 - x}^{\mu} - 1}x^{\nu} \over x}\,\dd x \,\right\vert_{\ \mu = 0\,,\ \nu\ =\ 0} \\[5mm] & = \left.\partiald[3]{}{\mu}\partiald{}{\nu} \int_{0}^{1}\bracks{\pars{1 - x}^{\mu}x^{\nu - 1} - x^{\nu - 1}}\dd x \,\right\vert_{\ \mu = 0\,,\ \nu\ =\ 0} \\[5mm] & = \partiald[3]{}{\mu}\partiald{}{\nu} \bracks{{\Gamma\pars{\mu + 1}\Gamma\pars{\nu} \over \Gamma\pars{\mu + \nu + 1}} - {1 \over \nu}}_{\ \mu = 0\,,\ \nu\ =\ 0} \\[5mm] & = \partiald[3]{}{\mu}\partiald{}{\nu} \braces{{\Gamma\pars{\mu + 1}\bracks{\Gamma\pars{\nu + 1}/\nu} \over \Gamma\pars{\mu + \nu + 1}} - {1 \over \nu}}_{\ \mu = 0\,,\ \nu\ =\ 0} \\[5mm] & = \partiald[3]{}{\mu}\partiald{}{\nu} \bracks{{\Gamma\pars{\mu + 1}\Gamma\pars{\nu + 1} - \Gamma\pars{\mu + \nu + 1} \over \nu\,\Gamma\pars{\mu + \nu + 1}}}_{\ \mu = 0\,,\ \nu\ =\ 0} \end{align} \begin{equation} \mbox{This limit is a 'cumbersome task'. Its value is}\ \bbx{\large\mc{I}_{2} = 12\,\zeta\pars{5} - \pi^{2}\zeta\pars{3}} \label{3}\tag{3} \end{equation}
\eqref{1}, \eqref{2} and \eqref{3} lead to
$$ \bbx{\sum_{k = 1}^{\infty}{H_{k} \over k^{4}} = 3\,\zeta\pars{5} - {1 \over 6}\,\pi^{2}\,\zeta\pars{3}} $$