Unfortunately I have no clever ideas to evaluate your integral. Instead I will follow your suggested line of attack and write $\tanh^{-1} (x) = \frac{1}{2} \ln \left (\frac{1 + x}{1 - x} \right ).$
Doing so the integral becomes \begin{align*} I &= \frac{1}{2} \int_0^1 \ln \left (\frac{1 + x}{1 - x} \right ) \frac{\ln x}{x(1 - x^2)} \, dx\\ &= \frac{1}{2} \int_0^1 \ln \left (\frac{1 - x^2}{(1 - x)^2} \right ) \frac{\ln x}{x(1 - x^2)} \, dx\\ &= \frac{1}{2} \int_0^1 \frac{\ln (1 - x^2) \ln x}{x(1 - x^2)} \, dx - \int_0^1 \frac{\ln (1 - x) \ln x}{x (1 - x^2)} \, dx\\ &= \frac{1}{2} I_1 - I_2. \end{align*}
The integral $I_1$
Enforcing the substitution $x \mapsto \sqrt{x}$ gives $I_1 = \frac{1}{4} \int_0^1 \frac{\ln (1 - x) \ln x}{x(1 - x)} \, dx.$ Making use of the generating function for the harmonic numbers $H_n$, namely $\sum_{n = 1}^\infty H_n x^n = - \frac{\ln (1 - x)}{1 - x},$ we have $I_1 = -\frac{1}{4} \sum_{n = 1}^\infty \int_0^1 x^{n - 1} \ln x \, dx.$
Noting that $\int_0^1 x^{n - 1} \ln x \, dx = -\frac{1}{n^2},$ a result that can be established using integration by parts, the integral becomes $I_1 = \frac{1}{4} \sum_{n = 1}^\infty \frac{H_n}{n^2}.$
The resulting sum, known as an Euler sum can be readily found (see here for example). Its value is $\sum_{n = 1}^\infty \frac{H_n}{n^2} = 2 \zeta (3).$ Thus $I_1 = \int_0^1 \frac{\ln (1 - x^2) \ln x}{x(1 - x^2)} \, dx = \frac{1}{2} \zeta (3).$
The integral $I_2$
From a partial fraction decomposition of $\frac{1}{x(1 - x^2)} = \frac{1}{x} + \frac{1}{2(1 - x)} + \frac{1}{2(1 + x)},$ the integral for $I_2$ can be written as \begin{align*} I_2 &= \int_0^1 \frac{\ln (1 - x) \ln x}{x} \, dx + \frac{1}{2} \int_0^1 \frac{x \ln x \ln (1 - x)}{1 - x} \, dx + \frac{1}{2} \int_0^1 \frac{x \ln x \ln (1 - x)}{1 + x} \, dx\\ &= I_\alpha + \frac{1}{2} I_\beta + \frac{1}{2} I_\gamma. \end{align*}
The first two integrals are relatively easy the find, the third is more difficult.
For the first \begin{align*} I_\alpha &= \int_0^1 \frac{\ln (1 - x) \ln x}{x} \, dx\\ &= -\text{Li}_2 (x) \ln x \Big{|}^1_0 + \int_0^1 \frac{\text{Li}_2 (x)}{x} \, dx\\ &= \int_0^1 \frac{\text{Li}_2 (x)}{x} \, dx\\ &= \text{Li}_3 (x) \Big{|}_0^1 = \text{Li}_3 (1) = \zeta (3). \end{align*}
For the second \begin{align*} I_\beta &= \int_0^1 \frac{x \ln x \ln (1- x)}{1- x} \, dx, \end{align*} we again make use of the generating function for the harmonic numbers. Doing so we have $I_\beta = - \sum_{n = 1}^\infty H_n \int_0^1 x^{n + 1} \ln x \, dx.$
By parts, it can be shown that $\int_0^1 x^{n + 1} \ln x \, dx = -\frac{1}{(n + 2)^2},$ thus $I_\beta = \sum_{n = 1}^\infty \frac{H_n}{(n + 2)^2}.$
To find the resulting sum we begin by shifting the summation index $n \mapsto n - 1$. Thus $I_\beta = \sum_{n = 2}^\infty \frac{H_{n - 1}}{(n + 1)^2}.$ Now from properties for the harmonic numbers $H_n = H_{n - 1} + \frac{1}{n}. \tag1$ Thus \begin{align*} I_\beta &= \sum_{n = 2}^\infty \frac{H_n}{(n + 1)^2} - \sum_{n = 2}^\infty \frac{1}{n(n + 1)^2}\\ &= \sum_{n = 1}^\infty \frac{H_n}{(n + 1)^2} - \sum_{n = 1}^\infty \frac{1}{n(n + 1)^2}. \end{align*}
In the first sum, shifting the summation index again by $n \mapsto n - 1$ and applying property (1) resulting in \begin{align*} \sum_{n = 1}^\infty \frac{H_n}{(n + 1)^2} &= \sum_{n = 2}^\infty \frac{H_n}{n^2} - \sum_{n = 2}^\infty \frac{1}{n^3} = \sum_{n = 1}^\infty \frac{H_n}{n^2} - \sum_{n = 1}^\infty \frac{1}{n^3} = 2 \zeta (3) - \zeta (3) = \zeta (3). \end{align*}
For the second sum, as $\frac{1}{n(n + 1)^2} = \frac{1}{n} - \frac{1}{n + 1} - \frac{1}{(n + 1)^2},$ we have \begin{align*} \sum_{n = 1}^\infty \frac{1}{n(n + 1)^2} &= \sum_{n = 1}^\infty \left (\frac{1}{n} - \frac{1}{n + 1} \right ) - \sum_{n = 1}^\infty \frac{1}{(n + 1)^2}. \end{align*} The first sum telescopes to $1$. For the second sum, it is $\sum_{n = 1}^\infty \frac{1}{(n + 1)^2} = \sum_{n = 2}^\infty \frac{1}{n^2} = \sum_{n = 1}^\infty \frac{1}{n^2} - 1 = \zeta (2) - 1.$ So $\sum_{n = 1}^\infty \frac{1}{n(n + 1)^2} = 1 - (\zeta (2) - 1) = 2 - \zeta (2),$ giving $I_\beta = \zeta (3) - (2 - \zeta (2)) = \zeta (3) + \zeta (2) - 2.$
The integral $I_\gamma$
Now for the difficult one. To evaluate it we will use double infinite sums. Since $\frac{\ln (1 - x)}{1 + x} = -\sum_{k = 0}^\infty \sum_{n = 0}^\infty \frac{(-1)^k}{n + 1} x ^{x + k + 1},$ the integral can be rewritten as $I_\gamma = - \sum_{k = 0}^\infty \sum_{n = 0}^\infty \frac{(-1)^k}{n + 1} \int_0^1 x^{n + k + 2} \ln x \, dx.$ As $\int_0^1 x^{n + k + 2} \ln x \, dx = - \frac{1}{(n + k + 3)^2},$ we can write $I_\gamma = \sum_{k = 0}^\infty \sum_{n = 1}^\infty \frac{(-1)^k}{(n + 1)(n + k + 3)^2}.$
By a partial fraction decomposition we have $\frac{1}{(n + 1)(n + k + 3)^2} = \frac{1}{(k + 2)^2(n + 1)} - \frac{1}{(k + 2)^2 (k + n + 3)} - \frac{1}{(k + 2)(k + n + 3)^2}.$ Now for the inner infinite sum over $n$ we have \begin{align*} \sum_{n = 0}^\infty \frac{1}{(n + 1)(n + k + 3)^2} &= \sum_{n = 0}^\infty \left [\frac{1}{(k + 2)^2(n + 1)} - \frac{1}{(k + 2)^2 (k + n + 3)} - \frac{1}{(k + 2)(k + n + 3)^2} \right ]\\ &= \frac{1}{(k + 2)^2} \sum_{n = 1}^{k + 2} \frac{1}{n} - \frac{1}{k + 2} \sum_{n = k + 3}^\infty \frac{1}{n^2}\\ &= \frac{1}{(k + 2)^2} H_{k + 2} - \frac{1}{k + 2} \sum_{n = 1}^\infty \frac{1}{n^2} + \frac{1}{k + 2} \sum_{n = 1}^{k + 2} \frac{1}{n^2}\\ &= \frac{H_{k + 2}}{(k + 2)^2} - \frac{\zeta (2)}{k + 2} + \frac{H^{(2)}_{k + 2}}{k + 2}. \end{align*} Here $H^{(a)}_n$ denotes the Generalised Harmonic numbers.
Thus $I_\gamma = \sum_{k = 0}^\infty \frac{(-1)^k H_{k + 2}}{(k + 2)^2} - \zeta (2) \sum_{k = 0}^\infty \frac{(-1)^k}{k + 2} + \sum_{k = 0}^\infty \frac{(-1)^k H^{(2)}_{k + 2}}{k + 2}.$
For the first sum, let $k \mapsto k - 2$, then $S_1 = \sum_{k = 2}^\infty \frac{(-1)^k H_k}{k^2} = 1 + \sum_{k = 1}^\infty \frac{(-1)^k H_k}{k^2} = 1 + A(1,2).$
Here $A(p,q) = \sum_{k = 1}^\infty \frac{(-1)^{k + 1} H^{(p)}_k}{k^q},$ correspond to the alternating Euler sums whose first few values can be found here. Since $A(1,2) = \frac{5}{8} \zeta (3)$ we have $S_1 = 1 - \frac{5}{8} \zeta (3).$
For the second sum $S_2 = \sum_{k = 0}^\infty \frac{(-1)^k}{k + 2} = -\sum_{k = 1}^\infty \frac{(-1)^{k + 1}}{k} + 1 = -\ln (2) + 1.$
For the third sum, let $k \mapsto k - 2$, then $S_3 = \sum_{k = 2}^\infty \frac{(-1)^k H^{(2)}_k}{k} = 1 - \sum_{k = 1}^\infty \frac{(-1)^{k + 1} H^{(2)}_k}{k} = 1 - A(2,1).$ And as $A(2,1) = \zeta (3) - \frac{1}{2} \zeta (2) \ln (2)$, we have $S_3 = 1 + \frac{1}{2} \zeta (2) \ln (2) - \zeta (3).$
So finally \begin{align*} I_\gamma &= \left (1 - \frac{5}{8} \zeta (3) \right ) - \zeta (2) \left (1 - \ln (2) \right ) + \left (1 + \frac{1}{2} \zeta (2) \ln (2) - \zeta (3) \right )\\ &= -\frac{13}{8} \zeta (3) + 2 - \zeta (2) + \frac{3}{2} \zeta (2) \ln (2). \end{align*}
So on combining all the results found we have \begin{align*} I_2 &= I_\alpha + \frac{1}{2} I_\beta + \frac{1}{2} I_\gamma\\ &= \zeta (3) + \frac{1}{2} \left (\zeta (3) + \zeta (2) - 2 \right ) + \frac{1}{2} \left (-\frac{13}{8} \zeta (3) + 2 - \zeta (2) + \frac{3}{2} \zeta (2) \ln (2) \right )\\ &= \frac{11}{16} \zeta (3) + \frac{3}{4} \zeta (2) \ln (2). \end{align*} And \begin{align*} I &= \frac{1}{2} I_1 - I_2\\ &= \frac{1}{4} \zeta (3) - \left (\frac{11}{16} \zeta (3) + \frac{3}{4} \zeta (2) \ln (2) \right )\\ &= -\frac{7}{16} \zeta (3) - \frac{3}{4} \zeta (2) \ln (2), \end{align*} or $\int_0^1 \frac{\tanh^{-1} (x) \ln (x)}{x(1 - x^2)} \, dx = -\frac{7}{16} \zeta (3) - \frac{\pi^2}{8} \ln (2),$ as required.