Evaluating $\sum_{n=1}^\infty\frac{(H_n)^2}{n}\frac{\binom{2n}n}{4^n}$ 
Question: How can we evaluate $$\sum_{n=1}^\infty\frac{(H_n)^2}{n}\frac{\binom{2n}n}{4^n},$$where $H_n=\frac11+\frac12+\cdots+\frac1n$?

Quick Results
This series converges because $$\frac{(H_n)^2}{n}\frac{\binom{2n}n}{4^n}=O\left(\frac{\ln^2n}{n^{3/2}}\right).$$
My Attempt
Recall the integral representation of harmonic number $$H_n=\int_0^1\frac{1-x^n}{1-x}d x$$
we have $$
S=\sum_{n=1}^\infty\frac1n\frac{\binom{2n}n}{4^n}\iint_{[0,1]^2}\frac{(1-x^n)(1-y^n)}{(1-x)(1-y)}d xd y\\
=\tiny\iint_{[0,1]^2}\frac{x y \log (4)-2 x y \log \left(\sqrt{1-x}+1\right)-2 x y \log \left(\sqrt{1-y}+1\right)+2 x y \log \left(\frac{1}{2} \left(\sqrt{1-x y}+1\right)\right)}{\left(\sqrt{1-x y}-1\right) \left(\sqrt{1-x y}+1\right)}dxdy\\
$$
This integral is too hard for me and Mathematica to compute. Numerical integration returns $12.6178$, it agrees with the numerical summation of the original series. I tried to integrate with respect to $x$, but failed.

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)$.
