How to evaluate, without contour integration the following integral: $$I=\int_0^1\frac{\ln^2x\ln(1+x)}{1+x^2}\ dx\ ?$$
@Cody mentioned in this solution that
$$I=\frac{\pi^{2}}{6}G+\frac{\pi^{3}}{32}\ln2-\frac{1}{768}\left[\psi_{3}\left(1/4\right)-\psi_{3}\left(3/4\right)\right]$$
but no proof was provided there, so any idea how to approach it?
The result from above can be further simplified by using digamma's reflection formula $$\psi(1-x)-\psi(x)=\pi\cot(\pi x)$$
And differentiating both sides three times then set $x=3/4$ to get
$$\psi_{3}(3/4)=16\pi^4-\psi_{3}(1/4)$$ $$\Rightarrow I=\frac{\pi^2}{6}G+\frac{\pi^{3}}{32}\ln2+\frac{\pi^4}{48}-\frac{1}{384}\psi_{3}(1/4)$$
Added:
Is it possible to evaluate $I$ using harmonic series?
$$\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)}$$