Evaluating $\int_0^1\ln(1+x^2)\ln(x^2+x^3)\frac{dx}{1+x^2}$

Let $a=\ln x, b=\ln(1-x), c=\ln(1+x), d=\ln(1+x^2)$. I use the following notations: $$I_{aa} = \int_0^1 \frac{\ln^2 x}{1+x^2}dx \qquad I_{ab} = \int_0^1 \frac{\ln x \ln(1-x)}{1+x^2}dx \qquad \cdots \qquad I_{cd} = \int_0^1 \frac{\ln (1+x) \ln(1+x^2)}{1+x^2}dx$$ Hence we get $10$ integrals. My goal is to find $9$ linearly independent relations between them, so your desired value $2I_{ad}+I_{cd}$ falls out easily.

Let $x=(1-u)/(1+u)$, then $dx/(1+x^2) = du/(1+u^2)$, and we have the following transformation rules: $$\begin{aligned}a &\mapsto b-c \\ b &\mapsto \ln 2 + a - c \\ c &\mapsto \ln 2 - c \\ d &\mapsto \ln 2 + d - 2c \end{aligned}$$

For example, we apply this on $I_{aa}$,we have $$\tag{1}I_{aa} = I_{bb} - 2I_{bc} + I_{cc}$$ We can apply this transformation on each of the ten integrals, but we only yield four linearly independent relations: $$\tag{2} I_{bb}=I_{aa}-2 I_{ac}-2 G \ln 2+I_{cc}$$ $$\tag{3} I_{dd}=2 \ln (2) \left(\frac{1}{2} \pi \ln (2)-G\right)+4 I_{cc}-4 I_{cd}+I_{dd}-\frac{1}{4} \pi \ln ^2(2)$$ $$\tag{4} I_{bd}=-2 I_{ac}+I_{ad}+\ln (2) \left(\frac{1}{2} \pi \ln (2)-G\right)-G \ln (2)+2 I_{cc}-I_{cd}-\frac{1}{8} \pi \ln ^2(2)$$

Of course, we have explicit evaluation of $I_{aa}$, which can be our fifth linearly independent relation: $$\tag{5} I_{aa} = \frac{\pi^3}{16}$$

To find more relations, we must rely on other methods. Here I use contour integration. Let $\log_1$ denote logarithm with branch cut at negative $x$-axis, while $\log_2$ denote logarithm with cut at positive $x$-axis. Integrate the function $$\frac{(\log_1 z)^a(\log_2 (z-1))^b}{1+z^2}$$ around a contour with two keyhole, wrapping around the two cuts: $(1,\infty)$ and $(-\infty,0)$. Then we obtain $$\int_1^\infty \cdots + \int_{-\infty}^0 \cdots = 2\pi i \text{(Sum of residues)}$$ The first integral's range can be brought back to $(0,1)$ via $x\mapsto 1/x$. The second integral, we first bring it back to $(0,\infty)$, then split intervals, finally apply $x\mapsto 1/x$ for the one with range $(1,\infty)$. After all these, We have $$\int_0^1 \frac{f_{a,b}(x)}{1+x^2} dx = 2\pi i \text{(Sum of residues)}$$ where $$f_{a,b}(x) = (-\ln (x))^a \left[(\ln (1-x)-\ln (x))^b-(\ln (1-x)-\ln (x)+2 \pi i)^b\right]-\left[(-\ln (x)-\pi i)^a-(-\ln (x)+\pi i)^a\right] (\ln (x+1)-\ln (x)+\pi i)^b-\left[(\ln (x)-\pi i)^a-(\ln (x)+\pi i)^a\right] (\ln (x+1)+\pi i)^b$$

Now apply this to $a=1,b=2$: $$\int_0^1 \frac{f_{1,2}(x)}{1+x^2}dx = -\frac{17 i \pi ^4}{16}+\frac{1}{4} i \pi ^2 \ln^2(2)-\pi ^3 \ln(2)$$ Hence comparing imaginary part:$$\tag{6}-2 \pi I_{aa}+4 \pi I_{ab}-4 \pi I_{ac}+4 \pi I_{cc}-\pi ^4=\frac{1}{4} \pi ^2 \ln ^2(2)-\frac{17 \pi ^4}{16}$$ This this our sixth linearly independent relation. Apply above method again to $a=0,b=3$: $$\tag{7}-6 \pi I_{bb}-6 \pi I_{aa}+12 \pi I_{ab}+2\pi^4 =-\frac{3}{4} \pi ^2 \ln (2)$$

The final two relations come from gamma/zeta function. Note that $$\int_1^\infty \frac{\ln^2(1+x^2)}{1+x^2}dx = I_{dd}-4I_{ad}+4I_{aa}$$ Hence $$\tag{8}2I_{dd}-4I_{ad}+4I_{aa} = \int_0^\infty \frac{\ln^2(1+x^2)}{1+x^2}dx = 4\int_0^{\pi/2} \ln^2(\cos x)dx = \frac{1}{6} \left(\pi ^3+12 \pi \ln ^2 2\right)$$

The last relation is more nontrivial: $$I_{ad}+I_{ab}+I_{ac} = \int_0^1 \frac{\ln x \ln \left(1-x^4 \right)}{1+x^2}dx = \frac{\pi^3}{16}-3G\ln 2 \tag{9}$$

which uses, in a critical way, values of digamma function.

Now solve those $9$ equations, we have one free variable (this involves a new constant, see below), and that free variable cancels for $2I_{ad}+I_{cd}$, proving your claim.

The new constant comes from $$\tag{10} I_{bb} = \int_0^1 \frac{\ln^2 x}{x^2-2x+2}dx = 2 \Im\left[\text{Li}_3\left(\frac{1+i}{2}\right)\right]$$

This follows directly from the indefinite integration: $$\int \frac{\ln^2 x}{x-a} = -2 \text{Li}_3\left(\frac{x}{a}\right)+2 \ln (x) \text{Li}_2\left(\frac{x}{a}\right)+\ln^2(x) \ln\left(1-\frac{x}{a}\right)$$

To consummate this approach, we obtain simultaneous evaluation of all $10$ integrals, all are nontrivial (except $I_{aa}, I_{bb}$) when considered individually. $$\begin{aligned} \int_0^1 \frac{\ln^2(1+x)}{1+x^2} dx &= -2 G \ln (2)-4 \Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)+\frac{7 \pi ^3}{64}+\frac{3}{16} \pi \ln ^2(2) \\ \int_0^1 \frac{\ln^2(1+x^2)}{1+x^2} dx &= -2 G \ln (2)+4 \Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)-\frac{7 \pi ^3}{96}+\frac{7}{8} \pi \ln ^2(2) \\ \int_0^1 \frac{\ln x \ln(1-x)}{1+x^2} dx &= \Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)-\frac{\pi ^3}{128}-\frac{1}{32} \pi \ln ^2(2) \\ \int_0^1 \frac{\ln x \ln(1+x)}{1+x^2} dx &= -2 G \ln (2)-3 \Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)+\frac{11 \pi ^3}{128}+\frac{3}{32} \pi \ln ^2(2) \\ \int_0^1 \frac{\ln x \ln(1+x^2)}{1+x^2} dx &= -G \ln (2)+2 \Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)-\frac{\pi ^3}{64}-\frac{1}{16} \pi \ln ^2(2) \\ \int_0^1 \frac{\ln (1-x) \ln(1+x)}{1+x^2} dx &= -G \ln (2)-\Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)+\frac{3 \pi ^3}{128}+\frac{3}{32} \pi \ln ^2(2) \\ \int_0^1 \frac{\ln (1-x) \ln(1+x^2)}{1+x^2} dx &= -\frac{1}{2} G \ln (2)+4 \Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)-\frac{5 \pi ^3}{64}+\frac{1}{8} \pi \ln ^2(2) \\ \int_0^1 \frac{\ln (1+x) \ln(1+x^2)}{1+x^2} dx &= -\frac{5}{2} G \ln (2)-4 \Im\left(\text{Li}_3\left(\frac{1+i}{2}\right)\right)+\frac{7 \pi ^3}{64}+\frac{3}{8} \pi \ln ^2(2) \end{aligned}$$

The Mathematica input is:

{aa -> \[Pi]^3/16, bb -> 2 Im[PolyLog[3, 1/2 + I/2]], cc -> (7 \[Pi]^3)/64 - 4 Im[PolyLog[3, 1/2 + I/2]] - 2 Catalan Log[2] + 3/16 \[Pi] Log[2]^2, dd -> -((7 \[Pi]^3)/96) + 4 Im[PolyLog[3, 1/2 + I/2]] - 2 Catalan Log[2] - 1/8 \[Pi] Log[2]^2 + 1/4 \[Pi] Log[4]^2, ab -> -(\[Pi]^3/128) + Im[PolyLog[3, 1/2 + I/2]] - 1/32 \[Pi] Log[2]^2, ac -> (11 \[Pi]^3)/128 - 3 Im[PolyLog[3, 1/2 + I/2]] - 2 Catalan Log[2] + 3/32 \[Pi] Log[2]^2, ad -> -(\[Pi]^3/64) + 2 Im[PolyLog[3, 1/2 + I/2]] - Catalan Log[2] - 1/16 \[Pi] Log[2]^2, bc -> (3 \[Pi]^3)/128 - Im[PolyLog[3, 1/2 + I/2]] - Catalan Log[2] + 3/32 \[Pi] Log[2]^2, bd -> -((5 \[Pi]^3)/64) + 4 Im[PolyLog[3, 1/2 + I/2]] - 1/2 Catalan Log[2] + 1/8 \[Pi] Log[2]^2, cd -> (7 \[Pi]^3)/64 - 4 Im[PolyLog[3, 1/2 + I/2]] - 5/2 Catalan Log[2] + 3/8 \[Pi] Log[2]^2}