Considering the algebraic identity
\begin{align*}
&(a-b)^3b = a^3b - 3a^2b^2 + 3ab^3 - b^4 = -2a^3b +3(a^3b+ab^3) -3a^2b^2 -b^4\\
&\Longrightarrow \ \ \ 2a^3b = -{b^4 \over 2} -{b^4 + 6a^2b^2\over 2} + 3(a^3b+ab^3) - (a-b)^3b
\end{align*} with $a = \ln(1-x)$ and $b= \ln (1+x)$ it follows that
\begin{align*}
2\int_0^1 {\ln^3(1-x)\ln(1+x)\over x}dx =& - \frac 1 2\int_0^1 {\ln^4(1+x)\over x}d x \\
&-\frac 12 \int_0^1 \frac{\ln^4(1+x) + 6\ln^2(1-x)\ln^2(1+x)}{x}dx\\
&+3\int_0^1 \frac{\ln^3(1-x)\ln(1+x) + \ln(1-x)\ln^3(1+x)}{x}dx\\
&- \int_0^1 \frac{\ln^3\left(\frac{1-x}{1+x}\right)\ln(1+x)}{x}dx\\
=:& -I_1 - I_2 + I_3 -I_4.
\end{align*}
For $I_1$, make substitution $y = \frac x {1+x}$ to get:
\begin{align*}
I_1 =& \frac 1 2 \int_0^{\frac 12} \frac{\ln^4(1-y)}{y(1-y)} dy \\
=& \frac 1 2\underbrace{ \int_0^{\frac 12} \frac{\ln^4(1-y)}{y} dy}_{z=1-y}+ \frac 1 2 \int_0^{\frac 12} \frac{\ln^4(1-y)}{1-y} dy\\
=& \frac 1 2 \int_{\frac 1 2 }^1 \frac{\ln^4 z} {1-z} dz + \frac {\ln^5 2}{10}\\
=& \frac 12 \sum_{n=1}^\infty \int_{\frac 1 2}^1 z^{n-1}\ln^4 z\ dz + \frac {\ln^5 2}{10}\\
=& \frac 12 \sum_{n=1}^\infty \frac{\partial^4}{\partial n^4}\left[\frac 1 n - \frac 1 {n2^n}\right] + \frac {\ln^5 2}{10}\\
=& \frac 12 \sum_{n=1}^\infty \left[\frac{24}{n^5} - \frac {24}{n^52^n} - \frac{24 \ln 2}{n^42^n}-\frac{12\ln^2 2}{n^3 2^n}-\frac{4\ln^3 2}{n^2 2^n} - \frac{\ln^4 2}{n2^n}\right] + \frac {\ln^5 2}{10}\\
=&12\zeta(5) - 12\text{Li}_5(1/2) - 12\ln 2 \text{Li}_4(1/2) -6\ln^2 2 \text{Li}_3(1/2) -2\ln^3 2\text{Li}_2(1/2)-\frac {2}{5}\ln^5 2\\
=&\boxed{-12\Big(\text{Li}_5(1/2) + \ln 2\text{Li}_4(1/2)-\zeta(5)\Big)-{21 \over 4}\zeta(3)\ln^2 2 +{1\over 3} \pi^2 \ln^3 2-{2 \over 5} \ln^5 2}
\end{align*} where the well-known values
\begin{align*}\text{Li}_2(1/2) = {\pi^2 \over 12}-{\ln^2 2\over 2} , \qquad \text{Li}_3(1/2) ={7\zeta(3) \over 8} -{\pi^2 \ln 2\over 12} + {\ln^3 2 \over 6}
\end{align*} are used.
Actually, $I_2$ was already evaluated by the OP here using the algebraic identity $$b^4 + 6a^2b^2 = \frac {(a-b)^4} 2+\frac{(a+b)^4}{2} -a^4.$$
It holds that
$$
\boxed{I_2 = \frac {21}{8} \zeta(5).}
$$
In fact, the value of $I_3$ can also be found in the previous answer of @Przemo's. For $I_3$, one can use the algebraic relation $3(a^3b + ab^3) =\frac 3 8 \left[ (a+b)^4 - (a-b)^4\right]$.
This gives
\begin{align*}
I_3=& \underbrace{\frac 3 8 \int_0^1 \frac{\ln^4(1-x^2)}{x} dx}_{x^2 = y} - \underbrace{\frac 3 8 \int_0^1 \frac{\ln^4\left(\frac{1-x}{1+x}\right)}{x} dx}_{\frac{1-x}{1+x} = y}\\
=&\frac 3 {16}\underbrace{\int_0^1 \frac{\ln^4(1-y)}{y} dy }_{1-y\mapsto y}- \frac 3 4 \int_0^1 \frac{\ln^4 y}{1-y^2} dy\\
=&\frac 3 {16}\int_0^1 \frac{\ln^4 y}{1-y} dy - \frac 3 4 \sum_{n=0}^\infty \int_0^1 y^{2n} \ln^4 y \ dy\\
=&\frac 3 {16}\sum_{n=1}^\infty \int_0^1 y^{n-1}\ln^4 y \ dy - \frac 3 4 \sum_{n=0}^\infty \frac {24}{(2n+1)^5}\\
=&\frac 3 {16}\sum_{n=1}^\infty \frac{24}{n^5} - 18 \sum_{n=0}^\infty \frac {1}{(2n+1)^5}\\
=&\frac {9}{2} \zeta(5)- 18\cdot \frac {31}{32}\zeta(5)\\
=&\boxed{-\frac{207}{16}\zeta(5)}
\end{align*} as can be found in @Przemo's answer.
For $I_4$, make substitution $ \frac{1-x}{1+x}\mapsto x$ to get
\begin{align*} I_4 = &2\int_0^1 \frac{\ln^3 x \ln\left(\frac 2 {1+x}\right)}{1-x^2} dx \\
=&2\ln 2 \int_0^1 \frac{\ln^3 x}{1-x^2} dx - \underbrace{2\int_0^1\frac{\ln^3 x \ln(1+x)}{1-x^2} dx }_{=:J}\\
=& 2\ln 2\sum_{n=0}^\infty \int_0^1 x^{2n} \ln^3 x\ dx - J\\
=& - 12\ln 2 \underbrace{\sum_{n=0}^\infty \frac 1 {(2n+1)^4}}_{\frac{15}{16}\zeta(4) = \frac{\pi^4}{96}} - J \\
=& -\frac{\pi^4 \ln 2}{8} - J.
\end{align*}
\begin{align*}
J = &\int_0^1\frac{2\ln^3 x \ln(1+x)}{1-x^2} dx \\
=& \underbrace{\int_0^1 \frac{\ln^3 x \ln(1+x)}{1+x}dx}_{=:A} + \int_0^1 \frac{\ln^3 x \ln(1+x)}{1-x}dx\\
=& A + \int_0^1 \frac{\ln^3 x \ln(1-x^2)}{1-x}dx -\int_0^1 \frac{\ln^3 x \ln(1-x)}{1-x}dx\\
=&A + \int_0^1 \frac{(1+x)\ln^3 x \ln(1-x^2)}{1-x^2}dx -\int_0^1 \frac{\ln^3 x \ln(1-x)}{1-x}dx\\
=&A + \underbrace{\int_0^1 \frac{\ln^3 x \ln(1-x^2)}{1-x^2}dx }_{=:B}+\underbrace{\int_0^1 \frac{x\ln^3 x \ln(1-x^2)}{1-x^2}dx}_{x^2 \mapsto x}-\int_0^1 \frac{\ln^3 x \ln(1-x)}{1-x}dx\\
=&A + B - \underbrace{\frac {15}{16} \int_0^1 \frac{\ln^3 x \ln(1-x)}{1-x}dx}_{=:C}\\
=&A + B - C.
\end{align*}
For $A$, we can use the McLaurin series of
$$
\frac{\ln (1+x)}{1+x} = \sum_{n=0}^\infty (-1)^{n-1}H_n x^n
$$ ($H_0= 0$) to get
\begin{align*}
A = & \sum_{n=0}^\infty (-1)^{n-1}H_n \int_0^1 x^n\ln^3 x \ dx \\
=&6 \sum_{n=0}^\infty \frac{(-1)^{n}H_n}{(n+1)^4}\\
=&6 \sum_{n=0}^\infty \frac{(-1)^{n}H_{n+1}}{(n+1)^4} - 6\sum_{n=0}^\infty \frac{(-1)^{n}}{(n+1)^5}\\
=&6 \sum_{n=1}^\infty \frac{(-1)^{n-1}H_{n}}{n^4} - 6\sum_{n=1}^\infty \frac{(-1)^{n-1}}{n^5}\\
=& 6\left(\frac{59}{32}\zeta(5) - \frac{\pi^2\zeta(3)}{12}\right)-6\cdot \frac{15}{16}\zeta(5)\\
=& \frac{87}{16}\zeta(5) - \frac{\pi^2 \zeta(3)}{2}.
\end{align*}
Here, the known value of $ \sum_{n=1}^\infty (-1)^{n-1}{H_n \over n^4}$ is used.
For $B$, make substitution $u = x^2$ to get
\begin{align*}
B =& \frac 1 {16} \int_0^1 \frac{\ln^3 u \ln(1-u)}{\sqrt u (1-u)} du \\
=& \frac 1 {16} \left[\frac{\partial^4}{\partial x^3\partial y} \text{B}(x,y)\right]_{x=\frac 1 2, y = 0^+}
\end{align*} where $\text{B}(\cdot,\cdot)$ is Euler's Beta function. We can use the fact that
\begin{align*}
\lim_{y\to 0^+}\frac{\partial^2}{\partial x\partial y} \text{B}(x,y) = -\frac 1 2 \psi''(x) + \psi'(x) \big[\psi(x) + \gamma\big]
\end{align*} to get
\begin{align*}
B =& \frac 1 {16}\frac{d^2}{dx^2}\left[-\frac 1 2 \psi''(x) + \psi'(x) \big[\psi(x) + \gamma\big]\right]_{x=\frac 1 2}\\
=&\frac 1 {16} \left[-\frac 1 2 \psi''''(1/2) + \psi'''(1/2)\big[\psi(1/2) + \gamma\big] + 3\psi'(1/2)\psi''(1/2)\right]\\
=& \frac 1 {16}\left[-21\pi^2 \zeta(3) + 372\zeta(5) - 2\pi^4 \ln 2\right]
\end{align*} which can be evaluated using the series representations of polygamma functions $$\psi(x) +\gamma = - \frac 1 x +\sum_{n=1}^\infty \frac 1 n - \frac 1 { n+x},\\
\psi'(x) = \sum_{n=0}^\infty \frac 1 {(n+x)^2}$$ and the derived fact that $\psi(\tfrac 1 2 )+\gamma = -2\ln 2$ and $\psi^{(k)}(\tfrac 1 2)=(-1)^{k+1}k!(2^{k+1}-1)\zeta(k+1)$ for $k\ge 1$.
For $C$, we can use the same method as used in the evaluation of $B$. It holds that
\begin{align*}
C =& \frac {15}{16} \left[\frac{\partial^4}{\partial x^3\partial y} \text{B}(x,y)\right]_{x=1, y = 0^+}\\
=&\frac {15} {16}\left[-\frac 1 2 \psi''''(1) + \psi'''(1)\big[\psi(1) + \gamma\big] + 3\psi'(1)\psi''(1)\right]\\
=&\frac{15}{16}\left[12\zeta(5) -6\zeta(2)\zeta(3)\right]\\
=&\frac {45}{4}\zeta(5) -\frac {15\pi^2 \zeta(3)}{16}
\end{align*} where $\psi(1) +\gamma = 0$, $\psi'(1) = \zeta(2)$, $\psi''(1) = -2\zeta(3)$ and $\psi''''(1) = -24\zeta(5)$ are used.
Combining $A,B,C$, we have that $$J =A+B-C= \frac{279}{16}\zeta(5) -\frac{7\pi^2\zeta(3)}{8} - \frac{\pi^4 \ln 2}{8}$$ and
$$
\boxed{I_4 = -\frac{\pi^4 \ln 2}{8} - J = -\frac{279}{16}\zeta(5)+\frac{7\pi^2\zeta(3)}{8}}
$$
Finally, these evaluate $\int_0^1 {\ln^3(1-x)\ln(1+x)\over x}dx =\frac 1 2\big[-I_1-I_2+I_3-I_4\big]$ as follows.
\begin{align*}
\int_0^1 {\ln^3(1-x)\ln(1+x)\over x}dx =&\ 6\text{Li}_5(1/2) + 6\ln 2\ \text{Li}_4(1/2)-\frac{81}{16}\zeta(5)-{7\pi^2 \over 16}\zeta(3)\\
&+\frac{21\ln^2 2}{8}\zeta(3)- \frac{1}{6}\pi^2\ln^3 2+\frac{1}{5}\ln^5 2.
\end{align*}
Using the identity given in the OP, we get the desired integral $I$
\begin{align*}
\int_0^{\frac 1 2}\frac{\text{Li}_2^2(x)}{x} dx = &-2\text{Li}_5(1/2) -2\ln 2\ \text{Li}_4(1/2)+\frac{27}{32}\zeta(5) +\frac{7\pi^2}{48}\zeta(3)-\frac{7\ln^2 2}{8}\zeta(3) \\
&-\frac{\pi^4\ln 2}{144} +\frac{\pi^2\ln^3 2}{12} - \frac{7\ln^5 2}{60}.
\end{align*}
I'll just show an idea that avoids those type of sum, but skip the calculations. You might also have better ideas to solve them.
For start we will denote $a=\ln(1-x)$ and $b=\ln(1+x)$ and use the following identity:
$$a^2=\frac12 (a+b)^2+\frac12(a-b)^2-b^2$$
$$\Rightarrow I=\frac12 \underbrace{\int_0^1 \frac{\ln x\ln^2(1-x^2)}{1+x}dx}_{I_1}+\frac12\underbrace{ \int_0^1 \frac{\ln x\ln^2\left(\frac{1-x}{1+x}\right)}{1+x}dx}_{I_2}-\underbrace{\int_0^1 \frac{\ln x\ln^2(1+x)}{1+x}dx}_{I_3}$$
For the first integral we will write the denominator as:
$$\frac{1}{1+x}=\frac{1}{1-x^2}-\frac{x}{1-x^2}$$
$$\Rightarrow I_1=\int_0^1 \frac{\ln x\ln^2(1-x^2)}{1-x^2}dx-{\int_0^1 \frac{x\ln x\ln^2(1-x^2)}{1-x^2}dx}$$
$$\overset{x^2\to x}=\frac14 \int_0^1 \frac{\ln x\ln^2(1-x)}{1-x}\frac{dx}{\sqrt x}-\frac14\int_0^1 \frac{\ln x\ln^2(1-x)}{1-x}dx$$
Those two integral can be found using the second identity from here.
Let's also take $I_2$ and substitute $\frac{1-x}{1+x}\to x$.
$$\Rightarrow I_2=\underbrace{\int_0^1 \frac{\ln(1-x)\ln^2 x}{1+x}dx}_{P}-\underbrace{\int_0^1 \frac{\ln(1+x)\ln^2 x}{1+x}dx}_{Q}$$
$$P-Q=I_2;\quad P+Q=\int_0^1 \frac{\ln(1-x^2)\ln^2 x}{1+x}dx$$
And again with the same trick done for $I_1$, we have:
$$P+Q=\int_0^1 \frac{\ln(1-x^2)\ln^2 x}{1-x^2}dx-\int_0^1 \frac{x\ln(1-x^2)\ln^2 x}{1-x^2}dx$$
$$\overset{x^2\to x}=\frac18\int_0^1 \frac{\ln(1-x)\ln^2 x}{1-x}\frac{dx}{\sqrt x}-\frac18 \int_0^1 \frac{\ln(1-x)\ln^2 x}{1-x}dx$$
Henceforth we can extract our second integral, $I_2$ as:
$$I_2=P-Q=(P+Q)-2Q$$
Note that $P+Q$ can again be found using the second identity from here.
Finally, we only need to find $Q$.
$$Q=\int_0^1 \frac{\ln(1+x)\ln^2 x}{1+x}dx=\sum_{n=1}^\infty (-1)^{n+1} H_n\int_0^1 x^{n}\ln^2 x=2\sum_{n=1}^\infty \frac{(-1)^{n+1}H_n}{(n+1)^3}$$
So $Q$ is actually an Euler sum in disguise, but you nicely found it here.
Also, $I_3$ is pretty easy, one just needs to use the same approach as done for $I_1$ in your following post.
$$I_3=\int_0^1 \frac{\ln x \ln^2(1+x)}{1+x}dx\overset{IBP}=-\frac12\int_0^1 \frac{\ln^3(1+x)}{x}dx$$
Best Answer
A solution in large steps by Cornel
By combining simple trigonometric formulae and Landen's Identity, we get
$$\mathcal{I}=\int_0^{\pi/2}\log(\sin(x))\operatorname{Li}_2(\sin^2(x))\textrm{d}x=\int_0^{\pi/2}\log(\sin(x))\operatorname{Li}_2\left(\frac{\tan^2(x)}{1+\tan^2(x)}\right)\textrm{d}x$$ $$=-\underbrace{\int_0^{\pi/2}\log(\sin (x))\operatorname{Li}_2(-\tan^2(x))\textrm{d}x}_{\mathcal{J}}-2\underbrace{\int_0^{\pi/2}\log(\sin(x))\log^2(\cos(x))\textrm{d}x}_{\text{Beta function: $\pi \zeta(3)/8-\pi\log ^3(2)/2$}} \tag1.$$
Now, based on the result in $(1.12)$ from the book, (Almost) Impossible Integrals, Sums, and Series, we have that $\displaystyle \operatorname{Li}_2(-\tan^2(x))=\int_0^1 \frac{\tan ^2(x) \log (y)}{1+\tan ^2(x)y} \textrm{d}y$, and then we write
$$\mathcal{J}=\int_0^{\pi/2}\left(\int_0^1 \frac{\log(\sin(x))\tan ^2(x) \log (y)}{1+\tan ^2(x)y} \textrm{d}y \right)\textrm{d}x=\int_0^1\left(\int_0^{\pi/2} \frac{\log(\sin(x))\tan ^2(x) \log (y)}{1+\tan ^2(x)y} \textrm{d}x \right)\textrm{d}y,$$ and if we make the change of variable $\tan(x)\mapsto x$, we arrive at $$\mathcal{J}=\int_0^1\left(\int_0^{\infty}\frac{x^2 \log (x)\log (y)}{(1+x^2)(1+y x^2)}\textrm{d}x\right)\textrm{d}y-2\int_0^1\left(\underbrace{\int_0^{\infty}\frac{y\log(y) x^2 \log(1+x^2) }{(1+x^2)(1+y^2 x^2)}\textrm{d}x}_{\displaystyle f(y)}\right)\textrm{d}y$$ or $$\mathcal{J}=\frac{\pi}{4}\int_0^1\frac{\log ^2(y)}{(y-1)\sqrt{y}}\textrm{d}y-2\log(2)\pi\int_0^1\frac{y\log (y)}{y^2-1}\textrm{d}y-2\pi\int_0^1\frac{\log^2 (y)}{y^2-1}\textrm{d}y$$ $$+\log(2)\pi\int_0^1\frac{\log (y)}{y-1}\textrm{d}y-\pi\int_0^1\frac{\log(1+y) \log (y)}{1+y}\textrm{d}y$$ $$+\pi\int_0^1\frac{\log \left(\frac{1+y}{2}\right) \log (y)}{y-1}\textrm{d}y$$ $$=\frac{\pi ^3}{6} \log (2)-\frac{7 }{8}\pi \zeta (3),\tag2$$ where all the resulting integrals are trivial, excepting the last one which is slightly more difficult, but it can be reduced to the calculations with Beta function, or one can use the generalization from A simple strategy of calculating two alternating harmonic series generalizations by Cornel Ioan Valean.
Combining $(1)$ and $(2)$ we arrive at the desired result.
A short note: For a fast calculation of $f(y)$ we may use the differentiation under the integral sign (with respect to $a$) e.g. $\displaystyle f(y,a)=\int_0^{\infty}\frac{y\log(y) x^2 \log(1+a^2 x^2) }{(1+x^2)(1+y^2 x^2)}\textrm{d}x$.
A BONUS: Using the integral $\mathcal{I}$, we observe that $$\mathcal{I}=\int_0^{\pi/2}\log(\sin(x))\operatorname{Li}_2(\sin^2(x))\textrm{d}x=\int_0^{\pi/2}\log(\sin(x))\operatorname{Li}_2(1-\cos^2(x))\textrm{d}x,$$ and if we combine the last integral with Dilogarithm reflection formula here, we get $$\mathcal{K}=\int_0^{\pi/2}\log(\sin(x))\operatorname{Li}_2(\cos^2(x))\textrm{d}x$$ $$=\zeta(2)\underbrace{\int_0^{\pi/2} \log (\sin(x))\textrm{d}x}_{-\log(2)\pi/2}-4\underbrace{\int_0^{\pi/2} \log ^2(\sin (x)) \log (\cos (x)) \textrm{d}x}_{\text{Beta function: $\pi \zeta(3)/8-\pi\log ^3(2)/2$}}$$ $$-\underbrace{\int_0^{\pi/2} \log (\sin (x))\operatorname{Li}_2(\sin ^2(x))\textrm{d}x}_{\displaystyle \mathcal{I}}$$ $$=\log ^3(2)\pi+\frac{1}{12}\log(2)\pi ^3-\frac{9 }{8}\pi \zeta (3). $$
The calculation of the integral $\mathcal{K}$ can also be achieved by similar steps to the ones used for the evaluation of the integral $\mathcal{I}$.
To conclude, we obtain that
Another short note: Alternatively, one can also build a system of relations with $\displaystyle \mathcal{I}=\int_0^{\pi/2}\log(\sin(x))\operatorname{Li}_2(\sin^2(x))\textrm{d}x$ and $\displaystyle \mathcal{K}=\int_0^{\pi/2}\log(\cos(x))\operatorname{Li}_2(\sin^2(x))\textrm{d}x$ and then calculate $\mathcal{I}+\mathcal{K}$ and $\mathcal{I}-\mathcal{K}$.