I will show that
$$I = \int_0^\infty \frac{x^2 \tan^{-1} x}{1 + x^2 + x^4} \, dx = \frac{\pi^2}{8 \sqrt{3}} + \frac{\pi}{24} \ln \left (\frac{2 - \sqrt{3}}{2 + \sqrt{3}} \right ) + \frac{2}{3} \mathbf{G},$$
where $\mathbf{G}$ is Catalan's constant.
Observing that
$$\frac{\tan^{-1} x}{x} = \int_0^1 \frac{dy}{1 + x^2 y^2},$$
on converting the integral to a double integral we have
$$I = \int_0^1 \int_0^\infty \frac{x^3}{(1 + x^2 + x^4)(1 + x^2 y^2)} \, dx \, dy,$$
after the order of integration has been changed. Finding a partial fraction decomposition for the integrand with respect to the variable $x$ gives
$$I = \int_0^1 \frac{1}{1 - y^2 + y^4} \int_0^\infty \left [- \frac{xy^2}{1 + x^2 y^2} + \frac{x - y^2 + 1}{2(x^2 + x + 1)} + \frac{x + y^2 - 1}{2(x^2 - x + 1)} \right ] \, dx \, dy.$$
Each of the $x$-integrals is elementary, being equal to either a log or the inverse tangent. The result is
$$I = \frac{\pi}{6 \sqrt{3}} \int_0^1 \frac{2y^2 - 1}{1 - y^2 + y^4} \, dy - \int_0^1 \frac{\ln y}{1 - y^2 + y^4} \, dy = \frac{\pi}{6\sqrt{3}} I_1 - I_2.$$
The first of the integrals is elementary. Here
\begin{align}
I_1 &= \int_0^1 \frac{2y^2 - 1}{1 - y^2 + y^4} \, dy\\
&= \int_0^1 \left [\frac{y - \sqrt{3}}{2 \sqrt{3} (-y^2 + y \sqrt{3} - 1)} + \frac{y + \sqrt{3}}{2 \sqrt{3} (y^2 + y \sqrt{3} + 1)} \right ] \, dy\\
&= -\frac{3}{4 \sqrt{3}} \int_0^1 \frac{2y + \sqrt{3}}{y^2 + y \sqrt{3} + 1} \, dy + \frac{1}{4} \int_0^1 \frac{dy}{(y + \sqrt{3}/2)^2 + 1/4}\\
& \qquad + \frac{3}{4 \sqrt{3}} \int_0^1 \frac{-2y + \sqrt{3}}{-y^2 + y \sqrt{3} - 1} \, dy + \frac{1}{4} \int_0^1 \frac{dy}{(y - \sqrt{3}/2)^2 + 1/4}\\
&=\frac{3}{4 \sqrt{3}} \ln \left (\frac{2 - \sqrt{3}}{2 + \sqrt{3}} \right ) + \frac{1}{2} \tan^{-1} (2 + \sqrt{3}) + \frac{1}{2} \tan^{-1} (2 - \sqrt{3})\\
&= \frac{\pi}{4} + \frac{3}{4 \sqrt{3}} \ln \left (\frac{2 - \sqrt{3}}{2 + \sqrt{3}} \right ).
\end{align}
The second of the integrals is a completely different beast. It will be handled through the dilogarithm machinery. In particular, we will make use of the result
$$\int_0^1 \frac{\ln x}{x - r} \, dx = \operatorname{Li}_2 \left (\frac{1}{r} \right ), \quad r \neq 0, \tag1$$
where $\operatorname{Li}_2 (x)$ is the dilogarithm. This result can be readily proved using integration by parts and noting the integral definition for $\operatorname{Li}_2 (x)$.
Denoting the roots of the equation $y^4 - y^2 + 1 = 0$ as $r_1 = (\sqrt{3} + i)/2$, $r_2 = -(\sqrt{3} + i)/2$, $r_3 = (\sqrt{3} - i)/2$, and $r_4 = (\sqrt{3} - i)/2$, on factoring the denominator appearing in the integrand of $I_2$ into linear factors over the complex domain before finding a partial fraction decomposition, one finds
\begin{align}
I_2 &= -\frac{\sqrt{3} + 3i}{12} \int_0^1 \frac{dx}{x - r_1} + \frac{\sqrt{3} + 3i}{12} \int_0^1 \frac{dx}{x - r_2}\\
& \qquad -\frac{\sqrt{3} - 3i}{12} \int_0^1 \frac{dx}{x - r_3} + \frac{\sqrt{3} - 3i}{12} \int_0^1 \frac{dx}{x - r_4},
\end{align}
or in terms of the dilogarithm, from (1)
\begin{align}
I_2 &= \frac{\sqrt{3} + 3i}{12} \left [\operatorname{Li}_2 \left (\frac{-\sqrt{3} + i}{2} \right ) - \operatorname{Li}_2 \left (\frac{\sqrt{3} - i}{2} \right ) \right ]\\
& \qquad + \frac{\sqrt{3} - 3i}{12} \left [\operatorname{Li}_2 \left (\frac{-\sqrt{3} - i}{2} \right ) - \operatorname{Li}_2 \left (\frac{\sqrt{3} + i}{2} \right ) \right ]\\
&= \frac{\sqrt{3} + 3i}{12} \left [\operatorname{Li}_2 \left (e^{5 \pi i/6} \right ) - \operatorname{Li}_2 \left (e^{-\pi i/6} \right ) \right ] + \frac{\sqrt{3} - 3i}{12} \left [\operatorname{Li}_2 \left (e^{-5 \pi i/6} \right ) - \operatorname{Li}_2 \left (e^{\pi i/6} \right ) \right ]\\
&= \frac{\sqrt{3} + 3i}{12} S_1 + \frac{\sqrt{3} - 3i}{12} S_2.
\end{align}
To find values for $S_1$ and $S_2$ the definition for the dilogarithm will be used, namely
$$\operatorname{Li}_2 (z) = \sum_{n = 1}^\infty \frac{z^n}{n^2}, \qquad |z| < 1.$$
For $S_1$, setting $z = e^{5\pi i/6}$ and $z = e^{-i\pi/6}$ in the above sum leads to
\begin{align}
S_1 &= \sum_{n = 1}^\infty \frac{1}{n^2} \left [\cos \left (\frac{5 n\pi}{6} \right ) - \cos \left (\frac{n\pi}{6} \right ) \right ] + i \sum_{n = 1}^\infty \frac{1}{n^2} \left [\sin \left (\frac{5 n\pi}{6} \right ) - \sin \left (\frac{n\pi}{6} \right ) \right ]\\
&= S_{1,1} + i S_{1,2}.
\end{align}
To find the first of these sums, as the series converges absolutely we can rearrange terms as follows:
\begin{align}
S_{1,1} &= - \sqrt{3} \sum_{\substack{n = 1\\n \in 1,13,\ldots}}^\infty \frac{1}{n^2} + \sqrt{3} \sum_{\substack{n = 1\\n \in 5,17,\ldots}}^\infty \frac{1}{n^2}\\
& \qquad + \sqrt{3} \sum_{\substack{n = 1\\n \in 7,19,\ldots}}^\infty \frac{1}{n^2} - \sqrt{3} \sum_{\substack{n = 1\\n \in 11,23,\ldots}}^\infty \frac{1}{n^2}\\
&= -\frac{\sqrt{3}}{144} \sum_{n = 0}^\infty \frac{1}{n + 1/12)^2} + \frac{\sqrt{3}}{144} \sum_{n = 0}^\infty \frac{1}{(n + 5/12)^2}\\
& \qquad + \frac{\sqrt{3}}{144} \sum_{n = 0}^\infty \frac{1}{(n + 7/12)^2} - \frac{\sqrt{3}}{144} \sum_{n = 0}^\infty \frac{1}{(n + 11/12)^2}\\
&= -\frac{\sqrt{3}}{144} \psi^{(1)} \left (\frac{1}{12} \right ) + \frac{\sqrt{3}}{144} \psi^{(1)} \left (\frac{5}{12} \right ) + \frac{\sqrt{3}}{144} \psi^{(1)} \left (\frac{7}{12} \right ) - \frac{\sqrt{3}}{144} \psi^{(1)} \left (\frac{11}{12} \right )\\
&= -\frac{\sqrt{3}}{144} \left [\psi^{(1)} \left (1 - \frac{1}{12} \right ) + \psi^{(1)} \left (\frac{1}{12} \right ) \right ]\\
& \qquad + \frac{\sqrt{3}}{144} \left [\psi^{(1)} \left (1 - \frac{5}{12} \right ) + \psi^{(1)} \left (\frac{5}{12} \right ) \right ]\\
&= -\frac{\sqrt{3}}{144} \frac{\pi^2}{\sin^2 (\pi/12)} + \frac{\sqrt{3}}{144} \frac{\pi^2}{\sin^2 (5\pi/12)}\\
&= - \frac{\pi^2}{6}.
\end{align}
Here $\psi^{(1)} (z)$ is the trigamma function with its reflexion formula having been used.
In a similar fashion the second sum $S_{1,2}$ can be found. Here
\begin{align}
S_{1,2} &= \sum_{\substack{n = 1\\n \in 1,13,\ldots}}^\infty \frac{1}{n^2} + 2 \sum_{\substack{n = 1\\n \in 3,15,\ldots}}^\infty \frac{1}{n^2}\\
& \qquad + \sum_{\substack{n = 1\\n \in 5,17,\ldots}}^\infty \frac{1}{n^2} - \sum_{\substack{n = 1\\n \in 7,19,\ldots}}^\infty \frac{1}{n^2}\\
& \qquad -2 \sum_{\substack{n = 1\\n \in 9,21,\ldots}}^\infty \frac{1}{n^2} - \sum_{\substack{n = 1\\n \in 11,23,\ldots}}^\infty \frac{1}{n^2}\\
&= \frac{1}{144} \sum_{n = 0}^\infty \frac{1}{(n + 1/12)^2} + \frac{1}{72} \sum_{n = 0}^\infty \frac{1}{(n + 3/12)^2}\\
& \qquad + \frac{1}{144} \sum_{n = 0}^\infty \frac{1}{(n + 5/12)^2} - \frac{1}{144} \sum_{n = 0}^\infty \frac{1}{(n + 7/12)^2}\\
& \qquad - \frac{1}{72} \sum_{n = 0}^\infty \frac{1}{(n + 9/12)^2} - \frac{1}{144} \sum_{n = 0}^\infty \frac{1}{(n + 11/12)^2}\\
&= \frac{1}{144} \left [\psi^{(1)} \left (\frac{1}{12} \right ) - \psi^{(1)} \left (\frac{7}{12} \right ) \right ] + \frac{1}{72} \left [\psi^{(1)} \left (\frac{1}{4} \right ) - \psi^{(1)} \left (\frac{3}{4} \right ) \right ]\\
& \qquad + \frac{1}{144} \left [\psi^{(1)} \left (\frac{5}{12} \right ) - \psi^{(1)} \left (\frac{11}{12} \right ) \right ]. \tag2
\end{align}
Two special values for the trigamma function are:
$$\psi^{(1)} \left (\frac{1}{4} \right ) = \pi^2 + 8 \mathbf{G} \quad \text{and} \quad \psi^{(1)} \left (\frac{3}{4} \right ) = \pi^2 - 8 \mathbf{G}.$$
The two other trigamma terms appearing in (2) can be dealt with by making use of the multiplication theorem for the trigamma function of
$$9 \psi^{(1)} (3z) = \psi^{(1)} (z) + \psi^{(1)} \left (z + \frac{1}{3} \right ) + \psi^{(1)} \left (z + \frac{2}{3} \right ).$$
Setting $z = 1/4$ leads to
\begin{align}
9 \psi^{(1)} \left (\frac{3}{4} \right ) &= \psi^{(1)} \left (\frac{1}{4} \right ) + \psi^{(1)} \left (\frac{7}{12} \right ) + \psi^{(1)} \left (\frac{11}{12} \right )\\
&= \psi^{(1)} \left (\frac{1}{4} \right ) + \psi^{(1)} \left (\frac{7}{12} \right ) + \psi^{(1)} \left (1 - \frac{1}{12} \right ),
\end{align}
or
$$\psi^{(1)} \left (\frac{1}{12} \right ) - \psi^{(1)} \left (\frac{7}{12} \right ) = 4 \sqrt{3} \pi^2 + 80 \mathbf{G},$$
where we have made use of the two special values for the trigamma function together with the reflexion formula.
Also, again starting with $z = 1/4$ in the multiplication formula we have
\begin{align}
9 \psi^{(1)} \left (\frac{3}{4} \right ) &= \psi^{(1)} \left (\frac{1}{4} \right ) + \psi^{(1)} \left (\frac{7}{12} \right ) + \psi^{(1)} \left (\frac{11}{12} \right )\\
&= \psi^{(1)} \left (\frac{1}{4} \right ) + \psi^{(1)} \left (1 - \frac{5}{12} \right ) + \psi^{(1)} \left (\frac{11}{12} \right ),
\end{align}
or
$$\psi^{(1)} \left (\frac{5}{12} \right ) - \psi^{(1)} \left (\frac{11}{12} \right ) = 4 \sqrt{3} \pi^2 + 80 \mathbf{G},$$
where again we have made use of the two special values for the trigamma function together with the reflexion formula. Thus
$$S_{1,2} = \frac{1}{144} \left (4 \sqrt{3} \pi^2 + 80 \mathbf{G} \right ) + \frac{1}{72} \left (\pi^2 + 8 \mathbf{G} -\pi^2 + 8 \mathbf{G}\right ) + \frac{1}{144} \left (-4 \sqrt{3} \pi^2 + 80 \mathbf{G} \right ) = \frac{4}{3} \mathbf{G}.$$
Thus
$$S_1 = -\frac{\pi^2}{6} + \frac{4 i \mathbf{G}}{3}.$$
Needless to say, a similar method can be used to find $S_2$. The result is:
$$S_2 = -\frac{\pi^2}{6} - \frac{4 i \mathbf{G}}{3},$$
giving
$$I_2 = \frac{\sqrt{3} + 3i}{12} \left (-\frac{\pi^2}{6} + \frac{4i \mathbf{G}}{3} \right ) + \frac{\sqrt{3} - 3i}{12} \left (-\frac{\pi^2}{6} - \frac{4i \mathbf{G}}{3} \right ) = -\frac{\pi^2}{12 \sqrt{3}} - \frac{2}{3} \mathbf{G}.$$
So on returning to our initial integral one has:
$$I = \frac{\pi}{6 \sqrt{3}} \left [\frac{\pi}{4} + \frac{3}{4 \sqrt{3}} \ln \left (\frac{2 - \sqrt{3}}{2 + \sqrt{3}} \right ) \right ] + \frac{\pi^2}{12 \sqrt{3}} + \frac{2}{3} \mathbf{G},$$
or
$$\int_0^\infty \frac{x^2 \tan^{-1} x}{1 + x^2 + x^4} \, dx = \frac{\pi^2}{8 \sqrt{3}} + \frac{\pi}{24} \ln \left (\frac{2 - \sqrt{3}}{2 + \sqrt{3}} \right ) + \frac{2}{3} \mathbf{G}.$$
This result agrees numerically with the previous result given by @Sangchul Lee in the comment section.
Best Answer
Write the integral as
$$I=\int\limits_0^{+\infty}\frac{e^{-ax^2}-e^{-bx^2}}{x}dx=\int\limits_0^{+\infty}\frac{e^{-(\sqrt{a}x)^2}-e^{-(\sqrt{b}x)^2}}{x}dx$$
Then, with $f(x)=e^{-x^2}$, apply the Frullani integral
$$I =\int\limits_0^{+\infty}\frac{f(\sqrt ax) -f(\sqrt bx)}{x}dx=[f(0)-f(\infty)]\ln \frac{\sqrt{b}}{\sqrt a}=\frac12\ln\frac ba$$