I was trying to compute the value of $$\sum_{k=1}^\infty \frac{H_k}{k^4}$$ and I was able to reduce it down to $$-\zeta(2)\zeta(3)+\int_0^1 \frac{\text{Li}_2^2(x)}{x}dx$$ However, I can't figure out how to compute the value of the integral $$\int_0^1 \frac{\text{Li}_2^2(x)}{x}dx$$ How can I find its value?
Don't try integration by parts. This integral is what I ended up with after integration by parts.
$\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}
\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}} $$
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)}. $$
Since $H_k=H_{k-1}+\frac1k$, we have $$\sum^\infty_{k=1}\frac{H_k}{k^4}=\sum^\infty_{k=1}\frac{H_{k-1}}{k^4}+\sum^\infty_{k=1}\frac1{k^5} =\sum^\infty_{k=1}\frac{H_{k-1}}{k^4}+\zeta(5).\tag1$$ If we define $$\sigma_h(s,t)=\sum^\infty_{n=1}\frac1{n^t}\sum^{n-1}_{k=1}\frac1{k^s},$$ already Euler expressed those sums by special values of the Riemann Zeta Function. He proved rigorously $$2\sigma_h(1,m)=m\,\zeta(m+1)-\sum^{m-2}_{k=1}\zeta(m-k)\,\zeta(k+1)$$ for $m\ge2$ (that's equation (3) in this paper). So the sum of the RHS of (1) is $\sigma_h(1,4)$, and the above formula easily gives $\sigma_h(1,4)=2\,\zeta(5)-\zeta(2)\,\zeta(3)$, and (1) gives the final result $$\sum^\infty_{k=1}\frac{H_k}{k^4}=3\,\zeta(5)-\zeta(2)\,\zeta(3).$$
To find your Euler sum I will make use of the following result which can be found here $$\sum_{k=1}^{\infty} \dfrac{H_k^{(p)}}{k^q}=\frac{(-1)^{q+1}}{\Gamma(q)} \int_{0}^{1} \frac{\ln^{q - 1} (u) \text{Li}_{p}(u)}{u (1 - u)}{du}.$$ Here $H^{(p)}_k$ are the generalised harmonic numbers with $H^{(1)}_k$ corresponding to the ordinary harmonic numbers $H_k$. It should be noted that the above result is proved using real methods only.
So when $p = 1$ and $q = 4$ we have $$\sum^\infty_{k = 1} \frac{H_k}{k^4} = -\frac{1}{\Gamma (4)} \int^1_0 \frac{\ln^3 (u) \, \text{Li}_1 (u)}{u (1 - u)} \, du = \frac{1}{6} \int^1_0 \frac{\ln^3 (u) \ln (1 - u)}{u (1 - u)} \, du,$$ where $\text{Li}_1 (u) = - \ln (1 - u)$ has been used.
The resulting integral can be found by reducing it to a double limit of the derivative of the beta function $\text{B}(x,y)$ as follows. As $$\text{B}(x,y) = \int^1_0 t^{x - 1} (1 - t)^{y - 1} \, dt,$$ we have $$\lim_{x \to 0^+} \lim_{y \to 0^+} \partial^3_x \partial_y \text{B}(x,y) = \int^1_0 \frac{\ln^3 (u) \ln (1 - u)}{u (1 - u)} \, du.$$ Thus \begin{align*} \sum^\infty_{k = 1} \frac{H_k}{k^4} &= \frac{1}{6} \lim_{x \to 0^+} \lim_{y \to 0^+} \partial^3_x \partial_y \text{B}(x,y) = \frac{1}{6} \left [-\frac{3}{4} \psi^{(4)} (1) + 3 \psi^{(2)}(1) \psi^{(1)}(1) \right ]. \end{align*} Here $\psi^{(m)}(x)$ is the polygamma function of order $m$. Values for the polygamma function when their arguments are equal to unity are well known. They are: $$\psi^{(4)}(1) = - 24 \zeta(5), \,\, \psi^{(2)}(1) = -2 \zeta (3), \,\,\psi^{(1)}(1) = \zeta (2),$$ and yields $$\sum^\infty_{k = 1} \frac{H_k}{k^4} = 3 \zeta (5) - \zeta (2) \zeta (3).$$
