Let the original integral equal $I$. Notice that
$$I = \int_{0}^{\infty}\frac{e^{4x}-e^{2x}}{x\left(e^{2x}+1\right)^{3}}dx = \int_{0}^{\infty}\frac{\tanh\left(x\right)}{4x\cosh^{2}\left(x\right)}dx$$
Integrating by parts yields
$$I = \frac{1}{8}\int_{0}^{\infty}\frac{\tanh^{2}\left(x\right)}{x^{2}}dx.$$
This popular link here provides a way to evaluate the integral above, but I'll add my own method.
Because the integrand is an even function, we can say
$$I = \frac{1}{16}\int_{-\infty}^{\infty}\frac{\tanh^{2}\left(x\right)}{x^{2}}dx.$$
Next, we will prove that
$$
\frac{1}{16}\int_{-\infty}^{\infty}\frac{\tanh^{2}\left(z\right)}{z^{2}}dz = \frac{1}{2\pi^{2}}\int_{0}^{1}\frac{1}{z}\ln^{2}\left(\frac{1-z}{1+z}\right)dz.
$$
Recall the Weierstrauss Product representation of $\cosh{(z)}$ as
$$
\cosh{(z)} = \prod_{n=0}^{\infty}\left(1+\left(\frac{2z}{\pi\left(2n+1\right)}\right)^{2}\right).
$$
Taking the logarithm and the derivative on both sides, we get
$$
\eqalign{
\ln\left(\cosh\left(z\right)\right)\ &=\ \ln\prod_{n=0}^{\infty}\left(1+\left(\frac{2z}{\pi\left(2n+1\right)}\right)^{2}\right) \\
&= \sum_{n=0}^{\infty}\ln\left(1+\left(\frac{2z}{\pi\left(2n+1\right)}\right)^{2}\right) \\
\frac{d}{dz}\ln\left(\cosh\left(z\right)\right)\ &=\ \frac{d}{dz}\sum_{n=0}^{\infty}\ln\left(1+\left(\frac{2z}{\pi\left(2n+1\right)}\right)^{2}\right) \\
\tanh\left(z\right) & = z\sum_{n=0}^{\infty}\frac{8}{4z^{2}+\left(\pi\left(2n+1\right)\right)^{2}} \\
\frac{\tanh\left(z\right)}{z}\ &=\ \sum_{n=0}^{\infty}\frac{8}{4z^{2}+\left(\pi\left(2n+1\right)\right)^{2}} \\
\frac{\tanh^{2}\left(z\right)}{z^{2}} &= \left(\sum_{n=0}^{\infty}\frac{8}{4z^{2}+\left(\pi\left(2n+1\right)\right)^{2}}\right)^{2} \\
&= \left(\sum_{n=0}^{\infty}\frac{8}{4z^{2}+\left(\pi\left(2n+1\right)\right)^{2}}\right)\left(\sum_{m=0}^{\infty}\frac{8}{4z^{2}+\left(\pi\left(2m+1\right)\right)^{2}}\right) \\
&= \sum_{n,m=0}^{\infty}\frac{64}{\left(4z^{2}+\left(\pi\left(2n+1\right)\right)^{2}\right)\left(4z^{2}+\left(\pi\left(2m+1\right)\right)^{2}\right)} \\
&= 64\sum_{n,m=0}^{\infty}\left(\frac{1}{4\left(z^{2}+\frac{\left(\pi\left(2n+1\right)\right)^{2}}{4}\right)}\frac{1}{4\left(z^{2}+\frac{\left(\pi\left(2m+1\right)\right)^{2}}{4}\right)}\right). \\
}
$$
On both sides, we integrate over $\mathbb{R}$ and multiply by $\frac{1}{16}$ so that $I$ becomes
$$
\eqalign{
&\frac{1}{16}\int_{-\infty}^{\infty}64\sum_{n,m=0}^{\infty}\left(\frac{1}{4\left(z^{2}+\frac{\left(\pi\left(2n+1\right)\right)^{2}}{4}\right)}\frac{1}{4\left(z^{2}+\frac{\left(\pi\left(2m+1\right)\right)^{2}}{4}\right)}\right)dz \\
&= \frac{2}{\pi^{2}}\sum_{n,m=0}^{\infty}\left(\frac{1}{\left(2n+1\right)\left(2m+1\right)}\frac{1}{2n+1+2m+1}\right) \\
&= \frac{2}{\pi^{2}}\sum_{n,m=0}^{\infty}\left(\frac{1}{\left(2n+1\right)\left(2m+1\right)}\int_{0}^{1}\frac{z^{2n+1+2m+1}}{z}dz\right) \\
&= \frac{2}{\pi^{2}}\int_{0}^{1}\frac{1}{z}\sum_{n,m=0}^{\infty}\frac{z^{2n+1}}{2n+1}\frac{z^{2m+1}}{2m+1}dz \\
&= \frac{1}{2\pi^{2}}\int_{0}^{1}\frac{1}{z}\ln^{2}\left(\frac{1-z}{1+z}\right)dz,\\
}
$$
where we used the power series representation
$$
\ln\left(\frac{1-z}{1+z}\right)=-2\sum_{k=0}^{\infty}\frac{z^{2k+1}}{2k+1},
$$
derived here in this YouTube video.
Next, we will prove that
$$
\int_{0}^{1}\frac{1}{z}\ln^{2}\left(\frac{1-z}{1+z}\right)dz = \frac{7}{2}\zeta{(3)},
$$
where Apery's Constant is defined as
$$
\zeta{(3)} = \sum_{n=1}^{\infty}\frac{1}{n^{3}}.
$$
Let that integral equal $J$ and let $f(z) = \frac{1}{z}\log^{2}\left(\frac{1-z}{1+z}\right)$. Notice we have the singularities $z_a = 0$ and $z_b = 1$. For small $\epsilon > 0$, let's define a closed contour $C$ as
$$
C = \left[1-\epsilon, \epsilon\right] \cup \gamma_1 \cup \left[i\epsilon, i\right] \cup \Gamma \cup \gamma_2.
$$
This contour travels clockwise. It closely resembles a quarter circle in the first quadrant of the complex plane such that $\gamma_1$ and $\gamma_2$ are small, quarter-circle, and indented arcs (with radius $\epsilon$) around the singularities $z_a = 0$ and $z_b = 1$, respectively. Also, $\Gamma$ closely resembles a quarter-circle arc (with radius $1$) such that a small radian, call it $t_{\epsilon}$, is cut off.
Below is a visual of $C$. One of my best friends Dairah coded this, so credit to her.

By Cauchy's Residue Theorem, we get $\oint_C f(z)dz$ to be
$$
\eqalign{
0 &= \int_{1-\epsilon}^{\epsilon}f(x)dx + \int_{\gamma_1}f(z)dz + \int_{i}^{i\epsilon}f(z)dz + \int_{\Gamma}f(z)dz + \int_{\gamma_2}f(z)dz \\
\int_{\epsilon}^{1-\epsilon}f(x)dx &= \int_{\gamma_1}f(z)dz + \int_{i}^{i\epsilon}f(z)dz + \int_{\Gamma}f(z)dz + \int_{\gamma_2}f(z)dz. \\
}
$$
Taking the limit as $\epsilon,t_{\epsilon} \to 0$ on both sides, we get
$$
J = \lim_{\epsilon \to 0}\int_{\gamma_1}f(z)dz + \lim_{\epsilon \to 0}\int_{i}^{i\epsilon}f(z)dz + \lim_{t_{\epsilon} \to 0}\int_{\Gamma}f(z)dz + \lim_{\epsilon \to 0}\int_{\gamma_2}f(z)dz.
$$
We will call each of those four limits $J_1$, $J_2$, $J_3$, and $J_4$, respectively.
First, we will prove $J_1 = 0$.
Let $z = \epsilon e^{i\phi}$ for $\phi \in \left[0, \frac{\pi}{2}\right]$. Then
$$
\eqalign{
J_1 &= \lim_{\epsilon \to 0}\int_{0}^{\frac{\pi}{2}} f\left(\epsilon e^{i\phi}\right)i \epsilon e^{i\phi}d\phi \\
&= i\lim_{\epsilon \to 0}\int_{0}^{\frac{\pi}{2}}\frac{1}{\epsilon e^{i\phi}}\log^2{\left(\frac{1-\epsilon e^{i\phi}}{1+\epsilon e^{i\phi}}\right)}\epsilon e^{i\phi}d\phi \\
&= i\int_{0}^{\frac{\pi}{2}}\log^2{\left(\frac{1}{1}\right)}d\phi \\
&= 0.
}
$$
Second, we will prove $J_4 = 0$ also.
Let $z = 1 + \epsilon e^{i\theta}$ where $\theta \in \left[\frac{\pi}{2}, \pi\right]$. Then using a chain of inequalities, observe that
$$
\eqalign{
\left|\int_{\frac{\pi}{2}}^{\pi}f\left(1+ \epsilon e^{i\theta}\right)i\epsilon e^{i\theta}\right| &= \left|\int_{\frac{\pi}{2}}^{\pi} \frac{1}{1+\epsilon e^{i\theta}}\log^2{\left(\frac{1-1-\epsilon e^{i\theta}}{1+1+\epsilon e^{i\theta}}\right)}i\epsilon e^{i\theta}d\theta\right| \\
&\leq \int_{\frac{\pi}{2}}^{\pi}\left|\frac{1}{1+\epsilon e^{i\theta}}\right|\cdot\epsilon\left|\log{\left(\frac{-\epsilon e^{i\theta}}{2+\epsilon e^{i\theta}}\right)}\right|^2 d\theta \\
&=\int_{\frac{\pi}{2}}^{\pi}\frac{\epsilon\left|\ln{\left|\frac{-\epsilon e^{i\theta}}{2+\epsilon e^{i\theta}}\right|}+i\arg{\left(\frac{-\epsilon e^{i\theta}}{2+\epsilon e^{i\theta}}\right)}\right|^2}{\sqrt{\epsilon^2+2\epsilon\cos{(\theta)}+1}} d\theta \\
&= \int_{\frac{\pi}{2}}^{\pi}\frac{\epsilon\left|\ln{\left(\frac{\epsilon}{\sqrt{\epsilon^2+4\epsilon\cos{(\theta)}+4}}\right)}+i\arg{\left(\frac{-\epsilon e^{i\theta}}{2+\epsilon e^{i\theta}}\right)}\right|^2}{\sqrt{\epsilon^2+2\epsilon\cos{(\theta)}+1}} d\theta \\
&= \int_{\frac{\pi}{2}}^{\pi}\frac{\left|\sqrt{\epsilon}\ln{(\epsilon)}-\sqrt{\epsilon}\ln{\left(\sqrt{\epsilon^2+4\epsilon\cos{(\theta)}+4}\right)}+\sqrt{\epsilon}i\arg{\left(\frac{-\epsilon e^{i\theta}}{2+\epsilon e^{i\theta}}\right)}\right|^2}{\sqrt{\epsilon^2+2\epsilon\cos{(\theta)}+1}} d\theta \\
}
$$
Taking $\epsilon \to 0$ of that last integral, we get
$$
\lim_{\epsilon \to 0}\int_{\frac{\pi}{2}}^{\pi}\frac{\left|\sqrt{\epsilon}\ln{(\epsilon)}-\sqrt{\epsilon}\ln{\left(\sqrt{\epsilon^2+4\epsilon\cos{(\theta)}+4}\right)}+\sqrt{\epsilon}i\arg{\left(\frac{-\epsilon e^{i\theta}}{2+\epsilon e^{i\theta}}\right)}\right|^2}{\sqrt{\epsilon^2+2\epsilon\cos{(\theta)}+1}} d\theta = \int_{\frac{\pi}{2}}^{\pi}\frac{\left|0-0+0\right|^2}{\sqrt{0+0+1}}d\theta = 0.
$$
(Notice that as $\epsilon \to 0$, we get $\sqrt{\epsilon}\ln{(\epsilon)} \to 0$ by L'Hôpital's Rule.)
So,
$$
0 \leq \lim_{\epsilon \to 0} \left|\int_{\frac{\pi}{2}}^{\pi}f\left(1+ \epsilon e^{i\theta}\right)i\epsilon e^{i\theta}\right| \leq 0.
$$
By the Squeeze Theorem, we get that
$$
\lim_{\epsilon \to 0} \left|\int_{\frac{\pi}{2}}^{\pi}f\left(1+ \epsilon e^{i\theta}\right)i\epsilon e^{i\theta}\right| = 0,
$$
implying that
$$
\lim_{\epsilon \to 0} \int_{\frac{\pi}{2}}^{\pi}f\left(1+ \epsilon e^{i\theta}\right)i\epsilon e^{i\theta} = 0.
$$
Third, we will prove $J_2 = \frac{7}{2}\zeta{(3)} - 2\pi C$, where $\zeta{(3)} := \sum_{n=0}^{\infty}\frac{1}{n^3}$ is Apery's Constant and $C := \sum_{n=0}^{\infty}\frac{(-1)^n}{(2n+1)^2}$ is Catalan's Constant.
We get $J_2$ to be
$$
\eqalign{
\int_{0}^{i}\frac{1}{z}\log^{2}\left(\frac{1-z}{1+z}\right)dz &= \int_{0}^{-1}\frac{1}{iz}\log^{2}\left(\frac{1-iz}{1+iz}\right)idz \text{ (using $z \to iz$)}\\
&= \int_{0}^{-1}\frac{\left(2i\arctan\left(z\right)\right)^{2}}{z}dz \\
&= 8\int_{0}^{1}\frac{\arctan\left(z\right)\ln\left(z\right)}{1+z^{2}}dz \text{ (using IBP)} \\
&= -4\pi\int_{1}^{\infty}\frac{\ln\left(z\right)}{z^{2}+1}dz+8\int_{1}^{\infty}\frac{\arctan\left(z\right)\ln\left(z\right)}{z^{2}+1}dz \text{ (using $z \to \frac{1}{z}$)} \\
&= -2\pi\int_{1}^{\infty}\frac{\ln\left(z\right)}{z^{2}+1}dz+4\int_{0}^{\infty}\frac{\arctan\left(z\right)\ln\left(z\right)}{z^{2}+1}dz, \\
}
$$
where we used the fact that $\int_0^1f(z)dz + \int_1^{\infty}f(z)dz = \int_0^{\infty}f(z)dz$ and did some algebra. Thus,
$$
-2\pi\int_{1}^{\infty}\frac{\ln\left(z\right)}{z^{2}+1}dz+4\int_{0}^{\infty}\frac{\arctan\left(z\right)\ln\left(z\right)}{z^{2}+1}dz = -2\pi C + 4\left(\frac{7}{8}\zeta(3)\right). \tag{1}
$$
Fourth, we will prove that $J_3 = 2\pi C$.
We parameterize $z = e^{it}$ for $t \in \left[t_{\epsilon}, \frac{\pi}{2}\right]$, for $t_{\epsilon} \to 0$. Then $J_3$ becomes
$$
\eqalign{
-\int_{0}^{\frac{\pi}{2}}\frac{1}{e^{it}}\log^{2}\left(\frac{1-e^{it}}{1+e^{it}}\right)ie^{it}dt &= -i\int_{0}^{\frac{\pi}{2}}\left(\ln\left|-i\tan\left(\frac{t}{2}\right)\right|+i\arg{\left(-i\tan\left(\frac{t}{2}\right)\right)}\right)^{2}dt \\
&= -2i\int_{0}^{\frac{\pi}{4}}\left(\ln\left(\tan\left(t\right)\right)-\frac{i\pi}{2}\right)^{2}dt \text{ (using $\frac{t}{2} \to t$)} \\
&= -2i\int_{0}^{\frac{\pi}{4}}\ln^{2}\left(\tan\left(t\right)\right)dt-2\pi\int_{0}^{\frac{\pi}{4}}\ln\left(\tan\left(t\right)\right)dt+\frac{i\pi^{3}}{8} \\
}
$$
Thus,
$$J_3 = -2i\left(\frac{\pi^{3}}{16}\right)-2\pi\left(-C\right)+\frac{i\pi^{3}}{8}. \tag{2} \\$$
Hence, $J$ simplifies to
$$
J = 0 -2\pi C + 4\left(\frac{7}{8}\zeta(3)\right) -2i\left(\frac{\pi^{3}}{16}\right)-2\pi\left(-C\right)+\frac{i\pi^{3}}{8} + 0 = \frac{7}{2}\zeta(3).
$$
Therefore, we multiply by $\frac{1}{2\pi^2}$ on both sides to solve for $I$.
In conclusion, the integral $I$ is
$$\int_{0}^{\infty}\frac{e^{4x}-e^{2x}}{x\left(e^{2x}+1\right)^{3}}dx = \frac{7\zeta(3)}{4\pi^2}.$$
If anyone wants to, I can solve $(1)$ and $(2)$ in the comments.
I greatly appreciate you taking the time to solve an integral I was trying to figure out for a long time. I attempted this integral a while back and used this attempt as part of another problem I was solving, so I've copied and pasted part of my solution from Overleaf and made a few adjustments. It is a long solution and there are plenty of easier ways to solve it like what Quanto provided in his answer, but I think using a crazy contour path is salvageable.
I hope I don't have complex argument/branch issues with my logarithms. I understand it is a lot to read, but please let me know if you catch any errors so I can learn more.