First, we prove a lemma on the integral representation of $(H_n)^2$.
$$I_n=\int_0^1\left(nx^{n-1}\ln^2(1-x)-\frac{x^n\ln x}{1-x}\right)d x-\zeta(2)=(H_n)^2$$
Let's prove by induction. $\displaystyle I_0=-\int_0^1\frac{\ln x}{1-x}dx=\zeta(2)=\zeta(2)+(H_0)^2$.\
Assume the equation holds for $n-1$,
$$\begin{aligned}
I_n&=\int_0^1\left(2(x^n-1)\frac{\ln(1-x)}{1-x}-\frac{x^n\ln x}{1-x}\right)d x-\zeta(2)\\
&=I_{n-1}+\int_0^1\left(2(x^n-x^{n-1})\frac{\ln(1-x)}{1-x}-\frac{(x^n-x^{n-1})\ln x}{1-x}\right)d x\\
&=(H_{n-1})^2+\int_0^1\left(-2x^{n-1}\ln(1-x)+x^{n-1}\ln x\right)d x\\
&=\left(H_n-\frac1n\right)^2-\frac1{n^2}+2\cdot\frac{H_n}n=(H_n)^2
\end{aligned}$$
Result Therefore, and by integrating $\displaystyle\sum_{n=1}^\infty\frac{\binom{2n}n}{4^n}x^n=\frac{1}{\sqrt{1-x}}-1$ from $0$ with respect to $x$, we have
$$\begin{aligned}
S&=\sum_{n=1}^\infty\frac1n\frac{\binom{2n}n}{4^n}\left(\int_0^1\left(nx^{n-1}\ln^2(1-x)-\frac{x^n\ln x}{1-x}\right)d x-\zeta(2)\right)\\
&=\int_0^1\left(\frac{1}{x\sqrt{1-x}}-\frac1x\right)\ln^2(1-x)d x-\int_0^12\ln\frac{2}{1+\sqrt{1-x}}\frac{\ln x}{1-x}d x-2\ln2\zeta(2)\\
&=I_1-I_2-2\ln2\zeta(2)
\end{aligned}$$
$I_1=12\zeta(3)$ can be easily deduced by substitution $x\mapsto 1-x^2$. $-2\ln2\zeta(2)+\frac32\zeta(3)$, the value of $I_2$, can also be deduced by the same substitution.
By combining these results, $S=\frac{21}2\zeta(3)$.
We proved in this solution that
$$\mathcal{I}=\int_0^1\frac{\tan^{-1}(x)\ln(1+x^2)}{x(1+x)}dx=\frac{\pi^3}{96}-\frac{\pi}{8}\ln^2(2)\tag1$$
and we proved here that
$$\tan^{-1}x\ln(1+x^2)=-2\sum_{n=1}^\infty\frac{(-1)^nH_{2n}}{2n+1}x^{2n+1}\tag2$$
By $(2)$ we get
$$\mathcal{I}=-2\sum_{n=1}^\infty\frac{(-1)^nH_{2n}}{2n+1}\int_0^1\frac{ x^{2n}}{1+x}dx$$
Using the identity $$\int_0^1\frac{x^{2n}}{1+x}dx=\ln2+H_n-H_{2n}$$
it follows that
$$\mathcal{I}=-2\ln(2)\sum_{n=1}^\infty\frac{(-1)^nH_{2n}}{2n+1}-2\sum_{n=1}^\infty\frac{(-1)^nH_{2n}H_n}{2n+1}+2\sum_{n=1}^\infty\frac{(-1)^nH_{2n}^{2}}{2n+1}$$
$$=-2\ln(2)\mathcal{S}_1-2\mathcal{S}_2+2\mathcal{S}_3\tag3$$
For $\mathcal{S}_1$ and $\mathcal{S}_3$, we use the classical identity:
$$\sum_{n=1}^\infty(-1)^n f(2n)=\Re\sum_{n=1}^\infty i^n f(n)$$
Therefore
$$\mathcal{S}_1=\Re\sum_{n=1}^\infty i^n\frac{H_n}{n+1}=\Re\left\{\frac{\ln^2(1-i)}{i}\right\}=-\frac{\pi}{8}\ln(2)\tag4$$
where we used $\sum_{n=1}^\infty x^n\frac{H_n}{n+1}=\frac{\ln^2(1-x)}{x}$ which follows from integrating $\sum_{n=1}^\infty H_nx^n=-\frac{\ln(1-x)}{1-x}$.
Similarly,
$$\mathcal{S}_3=\Re\sum_{n=1}^\infty i^n\frac{H_n^2}{n+1}$$
Using the generating function
$$\sum_{n=1}^\infty x^{n}\frac{ H_n^{2}}{n+1}=\frac{6\operatorname{Li}_3(1-x)-3\operatorname{Li}_2(1-x)\ln(1-x)-\ln^3(1-x)-3\zeta(2)\ln(1-x)-6\zeta(3)}{3x}$$
it follows that
$$\mathcal{S}_3=\Re\left\{\frac{6\operatorname{Li}_3(1-i)-3\operatorname{Li}_2(1-i)\ln(1-i)-\ln^3(1-i)-3\zeta(2)\ln(1-i)-6\zeta(3)}{3i}\right\}$$
$$=\Im\left\{\frac{6\operatorname{Li}_3(1-i)-3\operatorname{Li}_2(1-i)\ln(1-i)-\ln^3(1-i)-3\zeta(2)\ln(1-i)-6\zeta(3)}{3}\right\}$$
$$=2\Im\{\operatorname{Li}_3(1-i)\}+\frac12\ln(2)\ G+\frac{3\pi}{16}\ln^2(2)+\frac{5}{96}\pi^3\tag5$$
Plug the results of $(4)$ and $(5)$ in $(3)$ we get
$$\mathcal{I}=4\Im\{\operatorname{Li}_3(1-i)\}+\ln(2)\ G+\frac{5\pi}{8}\ln^2(2)+\frac{5}{48}\pi^3-2\mathcal{S}_2\tag6$$
By $(1)$ and $(6)$ we get
$$\mathcal{S}_2=\sum_{n=1}^\infty\frac{(-1)^nH_{2n}H_n}{2n+1}=2\Im\{\operatorname{Li}_3(1-i)\}+\frac12\ln(2)\ G+\frac{3\pi}{8}\ln^2(2)+\frac{3}{64}\pi^3$$
Best Answer
Let $\psi(z)$ be the digamma function.
The integral $$\int_{0}^{1} \frac{\ln(1+x^{2})}{1+x^{2}} \, \mathrm dx = -2 \int_{0}^{\pi/4} \log(\cos u) \, \mathrm du$$ can be evaluated in several ways. See this question for example.
But if you want to evaluate the series in a way that doesn't involve that integral, you could integrate the meromorphic function $$f(z) = \frac{\pi \csc(\pi z)\left(\psi(-z)+ \gamma \right)}{2z+1}$$ around a rectangular contour $C_{N}$ with vertices at $z = \pm (N+ \frac{1}{2}) \pm i \sqrt{N} $, where $N$ is a positive integer.
The function has double poles at the nonnegative integers, simple poles at the negative integers, and a simple pole at $z= - \frac{1}{2}$.
On the contour, the magnitude of $\psi(-z)$ grows slowly like $\ln|z|$ as $N \to \infty$.
This, combined with the fact that the magnitude of $\csc(\pi z)$ decays exponentially to zero as $\Im(z) \to \pm \infty$, means the integral will vanish on the contour as $N \to \infty$.
Therefore, we get $$ \lim_{N \to \infty} \int_{C_{N}} f(z) \, \mathrm dz = 0 = \sum_{n=-\infty}^{-1}\operatorname{Res}[f(z), n] + \sum_{n=0}^{\infty} \operatorname{Res}[f(z), n] + \operatorname{Res}[f(z), -1/2].$$
The Laurent series of $f(z)$ at a nonnegative integer is $$ \begin{align} f(z) &= \small \left(\frac{(-1)^n}{z-n} + \mathcal{O}(z-n) \right)\left(\frac{1}{z-n} + H_{n} + \mathcal{O}(z-n) \right) \left(\frac{1}{2n+1} -\frac{2(z-n)}{(2n+1)^{2}} + \mathcal{O}\left((z-n)^{2} \right)\right) \\ &= \frac{(-1)^{n}}{2n+1} \frac{1}{(z-n)^{2}} + \left( \color{red}{\frac{(-1)^{n}H_{n}}{2n+1} -\frac{2 (-1)^{n}}{(2n+1)^{2}} }\right ) \frac{1}{(z+n)}+ \mathcal{O}(1). \end{align}$$
Therefore, at an nonnegative integer $n$, there is a double pole with residue $$\frac{(-1)^{n}H_{n}}{2n+1} -\frac{2 (-1)^{n}}{(2n+1)^{2}}. $$
At $z=-n, n <0$, there is a simple pole with residue $$(-1)^{n+1} \frac{\psi(n)+ \gamma}{2n-1} = (-1)^{n+1} \frac{H_{n-1}}{2n-1}. $$
And at $z= - \frac{1}{2}$, there is a simple pole with residue $$-\frac{\pi}{2} \left( \psi \left(\frac{1}{2}\right)+ \gamma\right) = \pi \ln 2 . \tag{1}$$
Putting everything together, we have
$$\begin{align} 0 &= \sum_{n=1}^{\infty} \frac{(-1)^{n+1} H_{n-1}}{2n-1} + \sum_{n=0}^{\infty} \frac{(-1)^{n}H_{n}}{2n+1} -\sum_{n=0}^{\infty} \frac{2 (-1)^{n}}{(2n+1)^{2}} + \pi \ln 2 \\ &= -\sum_{k=0}^{\infty} \frac{(-1)^{k+1} H_{k}}{2k+1} - \sum_{n=0}^{\infty} \frac{(-1)^{n+1}H_{n}}{2n+1} -2 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)^{2}} + \pi \ln 2 \\ &= -2 \sum_{n=0}^{\infty}\frac{(-1)^{n+1}H_{n}}{2n+1} - 2 \sum_{n=0}^{\infty} \frac{ (-1)^{n}}{(2n+1)^{2}} + \pi \ln 2 \\ &= -2 \sum_{n=0}^{\infty}\frac{(-1)^{n+1}H_{n}}{2n+1} - 2 G + \pi \ln 2 \\ &= -2 \sum_{n=1}^{\infty}\frac{(-1)^{n+1}H_{n}}{2n+1} - 2 G + \pi \ln 2 \end{align}$$
The result then follows.
$(1)$ Special values $\psi \left(\frac12\right)$ and $\psi \left(\frac13\right)$