I tried 2 ways to find a closed form, though unsuccessful up to this point.
1st trial. Let $I$ denote the integral, and write
$$ I = 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \int_{0}^{\infty} \frac{e^{-x}(1 - e^{-(2n+1)x})}{x (1 + e^{-x})^{2}} \, dx. $$
In order to evaluate the integral inside the summation, we introduce new functions $I(s)$ and $J_n(s)$ by
$$ I(s) = 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \int_{0}^{\infty} \frac{x^{s-1} e^{-x}(1 - e^{-(2n+1)x})}{(1 + e^{-x})^{2}} \, dx =: 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} J_n(s) $$
so that $I = I(0)$. Then it is easy to calculate that for $\Re(s) > 1$, $J_n(s)$ is written as
$$ J_n(s) = \Gamma(s) \left( \eta(s-1) + \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s-1}} - (2n+1) \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s}} \right), $$
where $\eta$ is the Dirichlet eta function. Plugging this back to $I(s)$ and manipulating a little bit, we obtain
$$ I(s) = 8\Gamma(s) \left( \frac{\pi}{4} \eta(s-1) - 4^{-s}\left( \zeta(s, \tfrac{1}{4}) - \zeta(s, \tfrac{1}{2}) \right) + \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s-1}} \right). $$
This is valid for $\Re(s) > 1$. But if we can somehow manage to find an analytic continuation of the last summation part, then we may find the value of $I = I(0)$.
2nd trial. I began with the following representation
\begin{align*}
I
&= -2 \int_{0}^{\infty} \frac{1-e^{-t}}{1+e^{-t}} \left( \frac{1}{\cosh t} - \frac{2}{t} ( \arctan(1) - \arctan (e^{-t})) \right) \, \frac{dt}{t} \\
&= \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)(2n+2)} \int_{-\infty}^{\infty} \frac{\tanh^{2(n+1)} x}{x^{2}} \, dx.
\end{align*}
With some residue calculation, we can find that
\begin{align*}
\int_{-\infty}^{\infty} \frac{\tanh^{2n} x}{x^{2}} \, dx
&= \frac{2}{i\pi} \, \underset{z=0}{\mathrm{Res}} \left[ \psi_{1}\left(\tfrac{1}{2} + \tfrac{1}{i\pi} z\right) \coth^{2n} z \right] \\
&= 2^{2n+3} \sum_{m=1}^{n} (-1)^{m-1}m (1-2^{-2m-1}) A_{n-m}^{(2n)} \, \frac{\zeta(2m+1)}{\pi^{2m}},
\end{align*}
where $A_m^{(n)}$ is defined by the following combinatoric sum
$$ A_m^{(n)} = \sum_{\substack{ k_1 + \cdots + k_n = m \\ k_1, \cdots, k_n \geq 0 }} \frac{B_{2k_1} \cdots B_{2k_n}}{(2k_1)! \cdots (2k_n)!} = 2^{-2m} [z^{2m}](z \coth z)^{n} \in \Bbb{Q}, $$
where $B_k$ are Bernoulli numbers. Still the final output is egregiously complicated, so I stopped here.
3rd trial. The following yet another representation may be helpful, I guess.
$$ I = \int_{0}^{1/2} \frac{1 - \cot(\pi u/2)}{2} \left\{ \psi_1\left(\tfrac{1+u}{2}\right) - \psi_1\left(\tfrac{1-u}{2}\right) \right\} \, du. $$
Following the same approach as the given link, we have $\displaystyle\operatorname{arctanh}(x)=\frac12\ln\left(\frac{1+x}{1-x}\right)$, and with the same substitution $\displaystyle t=\frac{1+x}{1-x}$ we get
$$I(a,b)=-\frac1{2^{a-1}}\int_0^1\frac{(1+t)^{b-2}}{(1-t)^b}\ln^a(t)~\mathrm dt$$
Using generalized binomial expansion theorem, we can expand this as
$$I(a,b)=-\frac1{2^{a-1}(b-1)!}\sum_{j=0}^{b-2}\binom{b-2}j\sum_{k=0}^\infty(k+1)\dots(k+b-1)\int_0^1t^{k+j}\ln^a(t)~\mathrm dt$$
The integral at the end can be viewed as the $a$th derivative of $\displaystyle\int_0^1t^u~\mathrm dt=\frac1{u+1}$ at $u=k+j$, which is given by $\displaystyle\frac{(-1)^aa!}{(k+j+1)^{a+1}}$, which leaves us with the series
$$I(a,b)=\frac{(-1)^{a+1}a!}{2^{a-1}(b-1)!}\sum_{j=0}^{b-2}\binom{b-2}j\sum_{k=0}^\infty\frac{(k+1)\dots(k+b-1)}{(k+j+1)^{a+1}}$$
which boils the problem down to computing the last series. See that we have
$$\sum_{k=0}^\infty\frac{(k+1)\dots(k+b-1)}{(k+j+1)^{a+1}}=\sum_{k=j+1}^\infty\frac{P_{b,j}(k)}{k^{a+1}}=\sum_{m=0}^{b-1}\alpha_{b,j,m}\zeta(a-m+1)-\sum_{k=1}^j\frac{P_{b,j}(k)}{k^{a+1}}$$
for a polynomial $\displaystyle P_{b,j}(k)=\sum_{m=0}^{b-1}\alpha_{b,j,m}k^m$. This leaves us with the solution given by
$$I(a,b)=\frac{(-1)^{a+1}a!}{2^{a-1}(b-1)!}\sum_{j=0}^{b-2}\binom{b-2}j\left[\sum_{m=0}^{b-1}\alpha_{b,j,m}\zeta(a-m+1)-\sum_{k=1}^j\frac{P_{b,j}(k)}{k^{a+1}}\right]$$
The rest is just multiplying out $(k-j)\dots(k+b-j-2)$ to extract the coefficients of $P_{b,j}$. These are Stirling numbers of the first kind.
Best Answer
Integrate by parts to rewrite the integral as follows \begin{align} I_n=&\int_0^1 \frac{x^{n-1}\tanh^{-1} x}{1+x^n}dx\\ =& \frac1n \int_0^1 \tanh^{-1} x\ d\left(\ln\frac{1+x^n}2\right) \overset{ibp}=-\frac1n \int_0^1\frac{\ln\frac{1+x^n}2}{1-x^2}dx \end{align} Odd Case: Apply $$1+x^{2m+1}=(1+x)\prod_{k=1}^m (x^2+2x\cos a_k+1), \>\>\>\>\>a_k=\frac{2k\pi}{2m+1} $$ along with $\prod_{k=1}^m 2(1+\cos a_k)=1$ and $ \int_0^1\frac{\ln(1+t)}t dt=\frac{\pi^2}{12}$ \begin{align} I_{2m+1}=& \frac{-1}{2m+1} \int_0^1 \bigg(\ln\frac{1+x}2 + \sum_{k=1}^m \ln\frac{x^2+2x\cos a_k+1}{2+2\cos a_k}\bigg)\overset{x=\frac{1-t}{1+t}}{\frac{dx}{1-x^2}}\\ = & \ \frac{\pi^2}{24} -\frac1{2(2m+1)} \sum_{k=1}^m\int_0^1 \frac{\ln{(t^2 \tan^2\frac{a_k}2+1})}{t}dt \\ =& \ \frac{\pi^2}{24} +\frac1{4(2m+1)} \sum_{k=1}^m \text{Li}_2(-\tan^2\frac{a_k}2) \end{align}
Even Case: Analogous to the odd case $$ I_{2m} = \frac{\pi^2}{24} +\frac1{8m} \sum_{k=1}^{m} \text{Li}_2(-\tan^2\frac{b_k}2),\>\>\>\>\>b_k= \frac{(2k-1)\pi}{2m} $$ which can be further reduced per the symmetry of $b_k$, as well as the inverse formula $\text{Li}_2 (-z)+ \text{Li}_2 (-\frac1z)=-\frac{\pi^2}6-\frac12\ln^2 z $, i.e. \begin{align} I_{2m} = &\ \frac{\pi^2}{24}+\frac1{8m}\bigg[-\frac{m\pi^2}{12} -\frac12\sum_{k=1}^{[\frac m2]} \ln^2(\tan^2\frac{b_k}2)\bigg]\\ =&\ \frac{\pi^2}{32}-\frac1{4m}\sum_{k=1}^{[\frac m2]} \ln^2(\tan\frac{b_k}2) \end{align}