Closed form for the integral $\int_{0}^{\infty}\frac{\ln^{2}(x)\ln(1+x)}{(1-x)(x^{2}+1)}dx$
Integrate $f(z)=\dfrac{\ln^3{z}\ln(1+z)}{(1-z)(1+z^2)}$ along this contour.
The integral along the grey portions of the contour vanishes. We also need not worry about the removable singularity at $z=1$.
The imaginary part of the contour integral is \begin{align} {\rm Im}\oint_{\Gamma}f(z)\ {\rm d}z =&{\rm Im}\color{#E2062C}{\int^\infty_0\frac{\left(\ln^3{x}-(\ln{x}+2\pi i)^3\right)\ln(1+x)}{(1-x)(1+x^2)}{\rm d}x}+{\rm Im}\ \color{#6F00FF}{2\pi i\int^\infty_1\frac{\ln^3(-x)}{(1+x)(1+x^2)}{\rm d}x}\\ =&-6\pi\int^\infty_0\frac{\ln^2{x}\ln(1+x)}{(1-x)(1+x^2)}{\rm d}x+8\pi^3{\rm PV}\int^\infty_0\frac{\ln(1+x)}{(1-x)(1+x^2)}{\rm d}x\\ &+2\pi\int^\infty_1\frac{\ln^3{x}}{(1+x)(1+x^2)}{\rm d}x-6\pi^3\int^\infty_1\frac{\ln{x}}{(1+x)(1+x^2)}{\rm d}x \end{align} I will work out these integrals one by one. The first one is \begin{align} {\rm PV}\int^\infty_0\frac{\ln(1+x)}{(1-x)(1+x^2)}{\rm d}x =&\int^1_0\frac{\ln(1+x)}{1+x^2}{\rm d}x+\int^1_0\frac{x(1+x)\ln{x}}{1-x^4}{\rm d}x\\ =&\frac{\pi}{8}\ln{2}+\sum^\infty_{n=0}\int^1_0\left(x^{4n+1}+x^{4n+2}\right)\ln{x}\ {\rm d}x\\ =&\frac{\pi}{8}\ln{2}-\frac{1}{4}\sum^\infty_{n=0}\frac{1}{(2n+1)^2}-\sum^\infty_{n=0}\frac{1}{(4n+3)^2}\\ =&-\frac{\pi^2}{32}+\frac{\pi}{8}\ln{2}-\sum^\infty_{n=0}\frac{1}{(4n+3)^2} \end{align} Since \begin{align} \sum^\infty_{n=0}\frac{1}{(4n+1)^2}+\sum^\infty_{n=0}\frac{1}{(4n+3)^2}=&\frac{\pi^2}{8}\\ \sum^\infty_{n=0}\frac{1}{(4n+1)^2}-\sum^\infty_{n=0}\frac{1}{(4n+3)^2}=&\ \mathbf{G}\\ \end{align} The sum on the right evaluates to $\displaystyle \frac{\pi^2}{16}-\frac{\mathbf{G}}{2}$. So $${\rm PV}\int^\infty_0\frac{\ln(1+x)}{(1-x)(1+x^2)}{\rm d}x=\frac{\mathbf{G}}{2}-\frac{3\pi^2}{32}+\frac{\pi}{8}\ln{2}$$ The second one is \begin{align} \int^\infty_1\frac{\ln^3{x}}{(1+x)(1+x^2)}{\rm d}x =&-\int^1_0\frac{x(1-x)\ln^3{x}}{1-x^4}{\rm d}x\\ =&-\sum^\infty_{n=0}\int^1_0\left(x^{4n+1}-x^{4n+2}\right)\ln^3{x}\ {\rm d}x\\ =&\frac{3}{8}\sum^\infty_{n=0}\frac{1}{(2n+1)^4}-6\sum^\infty_{n=0}\frac{1}{(4n+3)^4}\\ =&\frac{\pi^4}{256}-\frac{1}{256}\psi_3\left(\frac{3}{4}\right) \end{align} The third one is \begin{align} \int^\infty_1\frac{\ln{x}}{(1+x)(1+x^2)}{\rm d}x =&-\int^1_0\frac{x(1-x)\ln{x}}{1-x^4}{\rm d}x\\ =&-\sum^\infty_{n=0}\int^1_0\left(x^{4n+1}-x^{4n+2}\right)\ln{x}\ {\rm d}x\\ =&\frac{1}{4}\sum^\infty_{n=0}\frac{1}{(2n+1)^2}-\sum^\infty_{n=0}\frac{1}{(4n+3)^2}\\ =&\frac{\mathbf{G}}{2}-\frac{\pi^2}{32} \end{align} Therefore, $${\rm Im}\oint_{\Gamma}f(z)\ {\rm d}z=-6\pi\int^\infty_0\frac{\ln^2{x}\ln(1+x)}{(1-x)(1+x^2)}{\rm d}x+\pi^3\mathbf{G}-\frac{71\pi^5}{128}-\frac{\pi}{128}\psi_3\left(\frac{3}{4}\right)+\pi^4\ln{2}$$ By the residue theorem, the contour integral is also equivalent to $2\pi i$ times the sum of residues. Keep in mind that $-\pi<\arg(1+z)\le\pi$ and $0\le\arg{z}<2\pi$. This means that $\ln(-i)=\small{\dfrac{i3\pi}{2}}$ but $\ln(1-i)=\small{\dfrac{1}{2}}\ln{2}-\small{\dfrac{i\pi}{4}}$. \begin{align} {\rm Im}\oint_{\Gamma}f(z)\ {\rm d}z =&{\rm Im}\ 2\pi i\left({\rm Res}(f,i)+{\rm Res}(f,-i)\right)\\ =&\frac{\pi^5}{64}-\frac{\pi^4}{32}\ln{2}-\frac{27\pi^5}{64}+\frac{27\pi^4}{32}\ln{2}\\ =&-\frac{13\pi^5}{32}+\frac{13\pi^4}{16}\ln{2} \end{align} Comparing both equalities, we get $$\int^\infty_0\frac{\ln^2{x}\ln(1+x)}{(1-x)(1+x^2)}{\rm d}x=\boxed{\displaystyle\Large{\color{#FF4F00}{\frac{\pi^2}{6}\mathbf{G}-\frac{1}{768}\psi_3\left(\frac{3}{4}\right)-\frac{19\pi^4}{768}+\frac{\pi^3}{32}\ln{2}}}}$$ which is equivalent to the proposed closed form.
Define $$ I=\int_{0}^{\infty}\frac{\ln^{2}(x)\ln(1+x)}{(1-x)(x^{2}+1)}dx, I(a)=\int_{0}^{\infty}\frac{\ln^{2}(x)\ln(1+ax)}{(1-x)(x^{2}+1)}dx, 0\le a\le 1. $$ Then $I(0)=0, I(1)=I$ and \begin{eqnarray} I'(a)&=&\int_{0}^{\infty}\frac{x\ln^{2}(x)}{(1+ax)(1-x)(x^{2}+1)}dx. \end{eqnarray} Define $$ f(z)=\frac{z}{(1+az)(1-z)(z^{2}+1)}. $$ Clearly $z=1$ is a removable singular point of $f(z)\ln^3 z$. Let $\Gamma$ to be the contour which is the line segment from $\varepsilon$ to $R$, together with two semicircles $S_\varepsilon$ and$S_R$ around 0 of radii $\varepsilon$, $R$ ($0<\varepsilon<1<R$). Clearly $f(z)$ is analytic inside $\Gamma$ except $z=\pm i,z=-\frac{1}{a}$ and $$ \text{Res}(f(z)\ln^3z,i)+\text{Res}(f(z)\ln^3z,-i)+\text{Res}(f(z)\ln^3z,-\frac{1}{a})=-\frac{\pi ^3 \left(a^2-1\right)+16 a (\pi i+\ln a)^3)}{16 \left(a^3+a^2+a+1\right)}. $$ It is easy to see $$ \bigg|\int_{S_\varepsilon}f(z)\ln^{3}(z)dz\bigg|\to 0 \text{ as }\varepsilon\to 0, \bigg|\int_{S_R}f(z)\ln^{3}(z)dz\bigg|\to 0 \text{ as }R\to\infty $$ and hence \begin{eqnarray} && \int_0^\infty f(x)\ln^{3}xdx-\int_0^\infty f(x)(\ln x+2\pi i)^3dx\\ &=&2\pi i(\text{Res}(f,i)+\text{Res}(f,-i)+\text{Res}(f(z)\ln^2z,\frac{1}{a}))\\ &=&-2\pi i\frac{\pi ^3 \left(a^2-1\right)+16 a (\pi i+\ln a)^3)}{16 \left(a^3+a^2+a+1\right)}. \end{eqnarray} Taking the imaginary parts for both sides gives $$ I'(a)=-\frac{16 a \ln a (\ln ^2a+\pi ^2)-3 \pi ^3 (a^2-1)}{48(a+1)(a^2+1)}. $$ Thus \begin{eqnarray} I(1)&=&-\int_0^1\frac{16 a \ln a (\ln ^2a+\pi ^2)-3 \pi ^3 (a^2-1)}{48(a+1)(a^2+1)}da\\ &=&-\int_0^1\frac{(1-a)[16 a \ln a(\ln ^2a+\pi ^2)-3 \pi ^3 (a^2-1)]}{48(1-a^4)}da\\ &=&-\int_0^1\frac{(1-a)a \ln a(\ln ^2a+\pi ^2)}{3(1-a^4)}da+3\int_0^1\frac{(1-a)\pi ^3 (a^2-1)}{16(1-a^4)}da\\ &=&-\frac{1}{3}\int_0^1\sum_{n=0}^\infty a^{4n}(1-a)a\ln a(\ln ^2a+\pi ^2)da+\frac{1}{64}\pi^3(-\pi+2\ln2)\\ &=&\sum_{n=0}^\infty\frac{1}{3}\left(\frac{6}{(4n+2)^4}-\frac{6}{(4n+3)^4}+\frac{\pi}{(4n+2)^2}-\frac{\pi}{(4n+3)^2}\right)+\frac{1}{64}\pi^3(-\pi+2\ln2)\\ &=&\frac{1}{768}(\psi_3(1/2)-\psi_3(3/4))+\frac{1}{48}(\psi_1(1/2)-\psi_1(3/4))+\frac{1}{64}\pi^3(-\pi+2\ln2). \end{eqnarray} It is well-know that \begin{eqnarray} \psi_3(1/2)=\pi^4, \psi_3(3/4)=8\pi^4-\beta(4),\psi_1(1/2)=\pi^2/2, \psi_1(3/4)=\pi^2-G, \end{eqnarray} and finally we have $$ I=\frac{G}{6}-\frac{1}{768}(19\pi^4-\psi_3(3/4))+\frac{1}{8}\pi^2\ln2. $$
I am going to go ahead and post my method. It is similar to xpauls except I used digamma, which is related to the harmonic series anyway.
Break integral up:
$$\int_{0}^{1}\frac{\log^{2}(x)\log(1+x)}{(1-x)(x^{2}+1)}dx+\int_{1}^{\infty}\frac{\log^{2}(x)\log(1+x)}{(1-x)(x^{2}+1)}dx$$
In the right integral, make the sub $x=1/t$. This gives:
$$\int_{0}^{1}\frac{\log^{2}(x)\log(1+x)}{(x^{2}+1)}dx+\int_{0}^{1}\frac{x\log^{3}(x)}{(1-x)(x^{2}+1)}dx$$
The right integral:
Break up into $$1/2\int_{0}^{1}\frac{x\log^{3}(x)}{x^{2}+1}dx-1/2\int_{0}^{1}\frac{\log^{3}(x)}{x^{2}+1}dx+1/2\int_{0}^{1}\frac{\log^{3}(x)}{1-x}dx$$
I am not going to work through each of these. But, suffice to say, they can be done without too much effort by using geometric series. For instance, take the middle one:
$$1/2\int_{0}^{1}\log^{3}(x)\sum_{k=0}^{\infty}(-1)^{k}x^{2k}dx=3\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k+1)^{4}}$$
Doing so to all three leads to series which evaluate in terms of $\zeta(4)$ and $\psi_{3}$. Summing them results in:
$$ \boxed{\displaystyle \int_{0}^{1}\frac{x\log^{3}(x)}{(1-x)(x^{2}+1)}dx=\frac{-9\pi^{4}}{256}+\frac{1}{512}\left[\psi_{3}(1/4)-\psi_{3}(3/4)\right]}$$
The left integral up top is a little more difficult. At least I think so.
$$\int_{0}^{1}\frac{\log^{2}(x)\log(1+x)}{x^{2}+1}dx$$
Use the Taylor series for $\log(1+x)$:
$$\int_{0}^{1}\frac{\log^{2}(x)}{x^{2}+1}\sum_{n=1}^{\infty}\frac{(-1)^{n}x^{n}}{n}$$
Note the incomplete Beta function defined as: $\displaystyle \int_{0}^{1}\frac{x^{a}}{x^{2}+1}dx=1/4\left[\psi \left(\frac{a+3}{4}\right)-\psi\left(\frac{a+1}{4}\right)\right]$.
Diffing this twice w.r.t 'a' introduces the log-square term and gives:
$$\int_{0}^{1}\frac{x^{a+n}\log^{2}(x)}{x^{2}+1}dx=1/64\left[\psi_{2} \left(\frac{a+n+3}{4} \right)-\psi_{2} \left(\frac{a+n+1}{4} \right) \right]$$.
Thus, letting $a=0$, $$\int_{0}^{1}\frac{\log^{2}(x)\log(1+x)}{x^{2}+1}dx=1/64\sum_{n=1}^{\infty}\frac{(-1)^{n}}{n}\left[\psi_{2}\left(\frac{n+3}{4}\right)-\psi_{2}\left(\frac{n+1}{4}\right)\right]$$
$$=\boxed{\displaystyle \frac{\pi^{2}}{6}G+\frac{\pi^{3}}{32}\log(2)-\frac{1}{768}\left[\psi_{3}\left(1/4\right)-\psi_{3}\left(3/4\right)\right]}$$
This series result, when combined with the other boxed result, gives the solution to the original integral.
The only minor issue I have is evaluating this tetragamma series. As I said, The Flajolet-Salvy residue method may work, but finding the correct kernel is the first important task. Since it alternates, I would assume something with $\pi \csc(\pi z)$
Of course, one could just say the heck with it and use this as a lemma. But, I would like to evaluate it though.