How to prove
$$\mathcal{I}=\int_0^1\frac{\tan^{-1}(x)\ln(1+x^2)}{x(1+x)}dx=\frac{\pi^3}{96}-\frac{\pi}{8}\ln^2(2)$$
This problem is proposed by a friend and I managed to compute it using only integration manipulation. My question is can we do it using harmonic series?
Here is my approach
Using the classical identity
$$\int_0^y\frac{\ln(1+yx)}{1+x^2}dx=\frac12 \tan^{-1}(y)\ln(1+y^2)$$
By integration by parts we have
$$\int_0^y\frac{y\tan^{-1}(x)}{1+yx}dx=\frac12 \tan^{-1}(y)\ln(1+y^2)$$
Now divide both sides by $y(1+y)$ then integrate from $y=0$ to $1$ we get
$$\frac12\mathcal{\color{red}{I}}=\frac12\int_0^1\frac{\tan^{-1}(y)\ln(1+y^2)}{y(1+y)}dy=\int_0^1\int_0^y\frac{\tan^{-1}(x)}{(1+y)(1+yx)}dxdy\\=\int_0^1\tan^{-1}(x)\left[\int_x^1\frac{dy}{(1+y)(1+yx)}\right]dx\\=\int_0^1\tan^{-1}(x)\left[\frac1{1-x}\ln\left(\frac{2(1+x^2)}{(1+x)^2}\right)\right]dx, \quad x=\frac{1-y}{1+y}\\=\int_0^1\left(\frac{\pi}{4}-\tan^{-1}(y)\right)\frac{\ln(1+y^2)}{y(1+y)}dy\\=\frac{\pi}{4}\int_0^1\frac{\ln(1+y^2)}{y(1+y)}dy-\mathcal{\color{red}{I}}$$
Rearranging the terms yields
$$\frac32\mathcal{\color{red}{I}}=\frac{\pi}{4}\int_0^1\frac{\ln(1+y^2)}{y(1+y)}dy=\frac{\pi}{4}\left(\frac{\pi^2}{16}-\frac34\ln^2(2)\right)\\ \Longrightarrow \mathcal{\color{red}{I}}=\frac{\pi^3}{96}-\frac{\pi}{8}\ln^2(2)$$
Where the last integral $\mathcal{J}=\int_0^1\frac{\ln(1+y^2)}{y(1+y)}dy$ can be proved using the common Feynman's method:
Let $$I(a)=\int_0^1\frac{\ln(1+a^2y^2)}{y(1+y)}dy, \quad I(0)=0, \quad I(1)=\mathcal{J}$$
$$I'(a)=\int_0^1\frac{2ay}{(1+y)(1+a^2y^2)}dy=2\frac{\tan^{-1}(a)}{1+a^2}+\frac{a\ln(1+a^2)}{1+a^2}-\ln(2)\frac{2a}{1+a^2}$$ $$\Longrightarrow \mathcal{J}=2\int_0^1\frac{\tan^{-1}(a)}{1+a^2}da+\int_0^1\frac{a\ln(1+a^2)}{1+a^2}da-\ln(2)\int_0^1\frac{2a}{1+a^2}da\\=\left(\frac{\pi}{4}\right)^2+\frac14\ln^2(2)-\ln^2(2)=\boxed{\frac{\pi^2}{16}-\frac34\ln^2(2)}$$