Here is what I have done
\begin{align} &\int_0^1\frac{\ln(x)\ln(1+x^2)}{1-x}dx\\ =&\int_0^1\frac{(1+x)(1+x^2) \log(x)\log(1-x^4)}{1-x^4} \ dx\\ &-\int_0^1\frac{(1+x)\log(x)\log(1-x^2)}{1-x^2} \ dx \end{align} and, then, after letting $x^4\mapsto x$ and $x^2\mapsto x$ respectively, use the beta function. Similarly we proceed with the second twin and we are done.
Please teach me another ways that only use real analysis.
Above also answers the first integral here
Denote $K_{\pm} = \int_0^1 \frac{\ln t}{1\pm t}dt$ and integrate the twins as follows \begin{align} I_\pm =&\int_0^1\frac{\ln x\ln(1+x^2)}{1\pm x}dx \\ =& \int_0^1 \ln(1+x^2)\> d\left(\int_0^x \frac{\ln t}{1\pm t}dt\right)\\ =& \>\ln2 \int_0^1 \frac{\ln t}{1\pm t}dt -\int_0^1 \frac{2x}{1+x^2}\left(\int_0^x \frac{\ln t}{1\pm t}\overset{t=xy}{dt}\right)dx\\ =& \>\ln2 K_{\pm}-2\int_0^1 \int_0^1\frac{x^2\ln(xy)}{(1+x^2)(1\pm xy)}dy\>dx\\ =& \>\ln2 K_{\pm} + 2\int_0^1 \int_0^1\frac1{1+y^2} \bigg( \frac{(1\mp xy)\ln(xy)}{1+x^2} -\frac{\overset{t=xy}{\ln(xy)}}{1\pm xy} \bigg)dx\>dy\\ =& \>\ln2 K_{\pm} \pm 2\int_0^1 \frac{x\ln x}{1+x^2} \overset{x^2\to x}{dx}\int_0^1\frac {2y}{1+y^2}dy\\ &+ 4\int_0^1 \frac{\ln x}{1+x^2} dx\int_0^1\frac 1{1+y^2}dy-\int_0^1 \frac{2}{y(1+y^2)}\left(\int_0^y \frac{\ln t}{1\pm t}dt\right) dy\\ = & \>\ln2 K_{\pm} \pm 2\left(\frac14K_+\right)\ln2 +4\>(-G)\>\frac\pi4 - \int_0^1 d\left(\ln\frac{y^2}{1+y^2}\right)\int_0^y \frac{\ln t}{1\pm t}dt\\ =& \>2\ln2 K_{\pm}\pm \frac12 \ln2 K_+ -\pi G+2\int_0^1 \frac{\ln^2y}{1\pm y}-I_{\pm}\\ \end{align} Substitute $K_- =-\frac{\pi^2}6$, $K_+ =-\frac{\pi^2}{12}$, $\int_0^1 \frac{\ln^2y}{1-y}dy=2\zeta(3)$, $\int_0^1 \frac{\ln^2y}{1+y}dy=\frac32\zeta(3)$ into above expression to arrive at \begin{align} &I_-=\int_0^1\frac{\ln x\ln(1+x^2)}{1-x}dx = -\frac{3\pi^2}{16}\ln2-\frac\pi2 G+2\zeta(3)\\ &I_+=\int_0^1\frac{\ln x\ln(1+x^2)}{1+x}dx =-\frac{\pi^2}{16}\ln2-\frac\pi2 G+\frac32\zeta(3) \end{align}