I am trying evaluating this $$J(k) = \int_{0}^{1} \frac{\ln^2x\ln\left ( \frac{1-x}{1+x} \right ) }{(x-1)^2-k^2(x+1)^2}\ \text{d}x.$$ For $k=1$, there has $$J(1)=\frac{\pi^4}{96}.$$ Maybe $J(k)$ doesn't have an explicit closed-form.
An integral relation: $$\int_{0}^{\infty} \frac{\arctan^3x}{1+k^2x^2} \text{d}x =-\frac{3}{2}J(k)+\frac{\pi^2}{8k} \left ( \operatorname{Li}_2\left ( \frac{1}{k} \right ) -\operatorname{Li}_2\left ( -\frac{1}{k} \right ) \right ) +\frac{\ln k}{8k} \ln^3\left ( \frac{k-1}{k+1} \right ),\qquad (k>1).$$
With some calculations, we followed that $$\int_{0}^{\infty} \frac{\arctan^3x}{1+k^2x^2} \text{d}x =\frac{\pi^4}{64k} +\frac{3}{4k}\left ( \operatorname{Li}_4\left ( \frac{k-1}{k+1} \right ) -\operatorname{Li}_4\left (- \frac{k-1}{k+1} \right ) \right ) +\frac{3\pi^2}{8k} \operatorname{Li}_2\left (- \frac{k-1}{k+1} \right ),\qquad(k>1).$$ Then the final result of $J(k)$ is $${\color{Green}{J(k) =\frac{\pi^2}{12k} \left ( \operatorname{Li}_2\left ( \frac{1}{k} \right ) -\operatorname{Li}_2\left ( -\frac{1}{k} \right ) \right ) +\frac{\ln k}{12k} \ln^3\left ( \frac{k-1}{k+1} \right ) -\frac{\pi^4}{96k} -\frac{1}{2k}\left ( \operatorname{Li}_4\left ( \frac{k-1}{k+1} \right ) -\operatorname{Li}_4\left (- \frac{k-1}{k+1} \right ) \right ) -\frac{\pi^2}{4k} \operatorname{Li}_2\left (- \frac{k-1}{k+1} \right ),\qquad (k>1).}}$$
I am not sure that a closed form is possible. However the problem can be made a bit nicer letting $$x=\frac{1-t}{1+t}\implies J(k)=-\frac 12\int_0^1 \frac {\log(t)}{k^2-t^2}\log ^2\left(\frac{1-t}{1+t}\right)\,dt$$ Now, using series expansion $$\frac{\log ^2\left(\frac{1-t}{1+t}\right)}{k^2-t^2}=\sum_{n=1}^\infty a_n\, k^{-2n}\,t^{2n}$$ with $$a_1=4 \qquad \text{and}\qquad a_n-a_{n-1}=4\,\frac{b_n}{c_n}\, k^{2n}$$ where the $b_n$ and $c_n$ correspond respectively to sequence $A002428$ and $A071968$ in $OEIS$ (use the absolute values for the $c_n$).
In fact, $\frac{b_n}{c_n}$ are the coefficients of the Taylor series of $\tanh ^{-1}(u)^2$.
Now, using integration by parts, $$\int_0^1 t^{2n} \, \log(t)\,dt=-\frac{1}{(2 n+1)^2}$$
From a programming point of view, the above seem to make things rather simple (except that the convergence is quite slow).
For example, for $k=2$ and using $1000$ terms in the summation, we have $J(2)\sim 0.092559$ while numerical integration lead to $J(2)=0.092561$.
Edit
Without any change of variable, we have, using the series expansion of $\log \left(\frac{1-x}{x+1}\right)$ and partial fraction decomposition
$$J(k)=\int_0^1\sum_{n=0}^\infty -\frac{2 \left(\frac{-k-1}{4 k ((k+1) x+k-1)}+\frac{k-1}{4 k ((k-1) x+k+1)}\right) x^{2 n+1} \log ^2(x)}{2 n+1}$$ which gives $$J(k)=\frac 1 k \sum_{n=0}^\infty \frac{\Phi \left(\frac{1-k}{1+k},3,2 n+1\right)-\Phi \left(\frac{1+k}{1-k},3,2 n+1\right)}{2 n+1}$$ where appears the Lerch transcendentfunction.
Now, we have a fast convergence. For $k=2$, the partial sums (up to $p$) are $$\left( \begin{array}{cc} p & \text{summation} \\ 0 & 0.0895113 \\ 1 & 0.0920160 \\ 2 & 0.0923839 \\ 3 & 0.0924835 \\ 4 & 0.0925206 \\ 5 & 0.0925374 \\ 6 & 0.0925460 \\ 7 & 0.0925509 \\ 8 & 0.0925539 \\ 9 & 0.0925558 \\ 10 & 0.0925570 \end{array} \right)$$ Using $a=\frac{1-k}{1+k}$, the summand $$b_n=\Phi \left(a,3,2 n+1\right)-\Phi \left(\frac1 a,3,2n+1\right)$$ and $$\Phi \left(a,3,2 n+1\right)=\frac 1{a^{2n+1}}\Bigg[\text{Li}_3(a)-\sum_{k=1}^{2n} \frac{a^k}{k^3} \Bigg]$$
Unfortunately, I was unable to reach the nice expression given by the OP in his/her edit.