Integrating $\int_0^1\frac{\ln^2x\ln(1+x)}{1+x^2} dx$ using real methods
$$\boxed{I=\int_0^1 \frac{\ln^2 x \ln(1+x)}{1+x^2}dx=\frac{\pi^3}{32}\ln 2 +\frac{\pi^2}{6}G-2\beta(4)}$$ $$\boxed{\int_0^1 \frac{\operatorname{Li}_3(-x)}{1+x^2}dx=\frac{\pi^2}{12} G-\frac{3\pi}{128}\zeta(3)-\beta(4)}$$ Tools used for the integral: $$\int_0^1 \frac{\ln x}{1+x^2}dx=-\beta(2)=-G\tag 1$$ $$\int_0^1 \frac{\ln^2 x}{1+x^2}dx=2\beta(3)=\frac{\pi^3}{16}\tag 2$$ $$ \int_0^1 \frac{\ln^3 x}{1+x^2}dx=-6\beta(4)=\frac{\pi^4}{16}-\frac{1}{128}\psi_3\left(\frac14\right)\tag 3$$ $$I'(a)=\int_0^\infty \frac{\ln^2 x}{(a+x)(1+x^2)}dx=-\frac{1}{3}\frac{\ln^3 a}{1+a^2}-\frac{\pi^2 }{3}\frac{\ln a}{1+a^2}+\frac{\pi^3}{8}\frac{a}{1+a^2}\tag 4$$ $\beta(x)$ is the Dirichlet Beta function and $G$ is Catalan's constant.
$(1)$,$(2)$ and $(3)$ follows easily by expanding the denominator into power series.
To prove $(4)$ we'll start by using partial fraction: $$I'(a)=\frac{1}{1+a^2}\int_0^\infty \ln^2 x\left(\frac{1}{a+x}-\frac{x}{1+x^2} \right)dx+\frac{a}{1+a^2}\int_0^\infty \frac{\ln^2 x}{1+x^2}dx$$ Now we will split the integrals in the point $1$, then use the substitution $x\to \frac{1}{x}$ in the $\int_1^\infty$ part and add it with the $\int_0^1$ part to get: $$I'(a)=\frac{1}{1+a^2}\int_0^1 \ln^2 x\left(\frac{1}{a+x}-\frac{1}{1/a+x}\right)dx+\frac{2a}{1+a^2}\int_0^1 \frac{\ln^2 x}{1+x^2}dx$$ $$=-\frac{2}{1+a^2}\left(\operatorname{Li}_3\left(-\frac{1}{a}\right)-\operatorname{Li}_3\left(-a\right)\right)+\frac{\pi^3}{8}\frac{a}{1+a^2}=-\frac{\ln a}{3}\frac{\pi^2 +\ln^2 a}{1+a^2}+\frac{\pi^3}{8}\frac{a}{1+a^2}$$ Above follows using one trilogarithm functional equation, which can be found here.
Evaluation of the integral: $$I=\int_0^1 \frac{\ln^2 x\ln(1+x)}{1+x^2}dx\overset{x\to \frac{1}{x}}=\int_1^\infty \frac{\ln^2 x(\ln(1+x)-\ln x)}{1+x^2}dx$$ $$\Rightarrow 2I=\int_0^\infty \frac{\ln^2 x \ln(1+x)}{1+x^2}dx-\int_1^\infty\frac{\ln^3 x}{1+x^2}dx$$
Now we will consider the following integral:
$$I(a)=\int_0^\infty \frac{\ln^2 x\ln(a+x)}{1+x^2}dx -\int_1^\infty \frac{\ln^3 x}{1+x^2}dx$$ Taking a derivative with respect to $a$ gives: $$ I'(a)= \int_0^\infty \frac{\ln^2 x}{(a+x)(1+x^2)}dx=-\frac{1}{3}\frac{\ln^3 a}{1+a^2}-\frac{\pi^2 }{3}\frac{\ln a}{1+a^2}+\frac{\pi^3}{8}\frac{a}{1+a^2}$$
We also have that:
$$I(0)=\int_0^\infty \frac{\ln^3 x}{1+x^2}dx-\int_1^\infty \frac{\ln^3 x}{1+x^2}dx=\int_0^1 \frac{\ln^3 x}{1+x^2}dx$$ And $2I= \left(I(1)-I(0)\right)+I(0)$, so: $$2I=-\frac13 \int_0^1 \frac{\ln^3 a}{1+a^2}da-\frac{\pi^2}{3}\int_0^1 \frac{\ln a}{1+a^2}da+\frac{\pi^3}{8}\int_0^1 \frac{a}{1+a^2}da+\int_0^1 \frac{\ln^3 x}{1+x^2}dx$$ $$=\frac{\pi^3}{16}\ln 2+\frac{\pi^2}{3}G -4\beta(4)\Rightarrow I=\boxed{\frac{\pi^3}{32}\ln 2 +\frac{\pi^2}{6}G-2\beta(4)}$$
Bonus: We can also obtain a nice result from this if we consider:
$$J(a)=\int_0^1 \frac{\ln^2 x\ln(1+ax)}{1+x^2}dx$$ $$\Rightarrow J'(a)=\frac{1}{1+a^2}\int_0^1\frac{x\ln^2 x}{1+x^2} dx+\frac{a}{1+a^2}\int_0^1 \frac{\ln^2 x}{1+x^2}dx-\frac{1}{1+a^2}\int_0^1\frac{\ln^2 x}{1/a+x}dx$$ $$=\frac{3\zeta(3)}{16}\frac{1}{1+a^2}+\frac{\pi^3}{16}\frac{a}{1+a^2}+\frac{2}{1+a^2}\operatorname{Li}_3(-a)$$ $$\Rightarrow I=\int_0^1 J'(a)da=\frac{3\pi}{64}\zeta(3)+\frac{\pi^3}{32}\ln 2+2\int_0^1 \frac{\operatorname{Li}_3(-a)}{1+a^2}da$$ And from this we can extract the last integral: $$\boxed{\int_0^1 \frac{\operatorname{Li}_3(-x)}{1+x^2}dx=\frac{\pi^2}{12} G-\frac{3\pi}{128}\zeta(3)-\beta(4)}$$
Let $$I=\int_0^1\frac{\ln^2x\ln(1+x)}{1+x^2}\ dx$$
and $$I(a)=\int_0^1\frac{\ln^2x\ln(1+ax)}{1+x^2}\ dx$$
$$I(0)=0 , \quad I(1)=I$$
\begin{align} I^{'}(a)&=\int_0^1\frac{x\ln^2x}{(1+ax)(1+x^2)}\ dx=\int_0^1\frac{\ln^2x}{1+a^2}\left(\frac{a}{1+x^2}+\frac{x}{1+x^2}-\frac{a}{1+ax}\right)\ dx\\ &=\frac{a}{1+a^2}\int_0^1\frac{\ln^2x}{1+x^2}\ dx+\frac1{1+a^2}\int_0^1\frac{x\ln^2x}{1+x^2}\ dx-\frac{a}{1+a^2}\int_0^1\frac{\ln^2x}{1+ax}\ dx\\ &=\frac{a}{1+a^2}\left(\frac{\pi^3}{16}\right)+\frac1{1+a^2}\left(\frac{3}{16}\zeta(3)\right)-\frac{a}{1+a^2}\left(\frac{-2\operatorname{Li}_3(-a)}{a}\right) \end{align}
Then
\begin{align} I&=\frac{\pi^3}{16}\int_0^1\frac{a\ da}{1+a^2}+\frac{3}{16}\zeta(3)\int_0^1\frac{da}{1+a^2}+2\int_0^1\frac{\operatorname{Li}_3(-a)}{1+a^2}\ da\\ &=\frac{\pi^3}{32}\ln 2+\frac{3\pi}{64}\zeta(3)+2\color{blue}{\int_0^1\frac{\operatorname{Li}_3(-a)}{1+a^2}\ da}\tag{1}\\ \end{align}
.
To evaluate the blue integral, I am going to use the same approach in my post here:
\begin{align} \int_0^1 \frac{\operatorname{Li}_3(-a)}{1+a^2}\ da&=\int_0^\infty \frac{\operatorname{Li}_3(-a)}{1+a^2}\ da-\underbrace{\int_1^\infty \frac{\operatorname{Li}_3(-a)}{1+a^2}\ da}_{a\mapsto 1/a}\\ &=\int_0^\infty \frac{\operatorname{Li}_3(-a)}{1+a^2}\ da-\int_0^1 \frac{\operatorname{Li}_3(-1/a)}{1+a^2}\ da\\ &\left\{\text{add the integral to both sides}\right\}\\ 2\int_0^1 \frac{\operatorname{Li}_3(-a)}{1+a^2}\ da&=\int_0^\infty\frac{\operatorname{Li}_3(-a)}{1+a^2}\ da-\int_0^1 \frac{\operatorname{Li}_3(-1/a)-\operatorname{Li}_3(-a)}{1+a^2}\ da\\ &\{\text{use}\ \operatorname{Li}_3(-1/a)-\operatorname{Li}_3(-a)=\frac16\ln^3a+\zeta(2)\ln a\}\\ &=\int_0^\infty\frac{\operatorname{Li}_3(-a)}{1+a^2}\ da-\frac16\underbrace{\int_0^1\frac{\ln^3a}{1+a^2}\ da}_{-6\beta(4)}-\zeta(2)\underbrace{\int_0^1\frac{\ln a}{1+a^2}\ da}_{-G}\\ &=\color{red}{\int_0^\infty\frac{\operatorname{Li}_3(-a)}{1+a^2}\ da}+\beta(4)+\zeta(2)G\tag{2} \end{align}
Also the red integral can be evaluated the same way in this solution:
$$\operatorname{Li}_{3}(-a)=\frac12\int_0^1\frac{-a\ln^2 u}{1+au}\ du$$
\begin{align} \int_0^\infty\frac{\operatorname{Li}_{3}(-a)}{1+a^2}\ da&=-\frac12\int_0^1\ln^2u\left(\int_0^\infty\frac{a}{(1+ua)(1+a^2)}\ da\right)\ du\\ &=-\frac12\int_0^1\ln^2u\left(\frac12\left(\frac{\pi u}{1+u^2}-\frac{2\ln u}{1+u^2}\right)\right)\ du\\ &=-\frac{\pi}{4}\underbrace{\int_0^1\frac{u\ln^2u}{1+u^2}\ du}_{\frac{3}{16}\zeta(3)}+\frac12\underbrace{\int_0^1\frac{\ln^3u}{1+u^2}\ du}_{-6\beta(4)}\\ &=\frac{-3\pi}{64}\zeta(3)-3\beta(4)\tag{3} \end{align}
Plugging (3) in (2) we get
$$2\int_0^1 \frac{\operatorname{Li}_3(-a)}{1+a^2}\ da=\zeta(2)G-\frac{3\pi}{64}\zeta(3)-2\beta(4)\tag{4}$$
Now plug (4) in (1) we get
$$I=\int_0^1\frac{\ln^2x\ln(1+x)}{1+x^2}\ dx=\frac{\pi^3}{32}\ln2+\zeta(2)G-2\beta(4)$$