How to calculate $\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx$? Using the derivative of beta function, find 
$$I=\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx$$
setting $x^2=y$ gives
$$I=\frac18\int_0^1\frac{\ln^2y\ln^2(1-y)}{\sqrt{y}(1-y)}\ dy=\frac18\left.\frac{\partial^4}{\partial a^2\partial b^2}\text{B}(a,b)\right|_{a\mapsto 1/2\\b\mapsto0^{+}}$$
Any good software that can find the 4th derivative and also gives the final result? Wolfram fails to calculate it (or maybe I do not know how to use it well) and when I tried to do it manually, some terms involve $\psi(b)$ and if we take the limit, then $\psi(0)$ is undefined and even if I take the limit of $\psi(b)$ together with other terms, still undefined. I do not know how to avoid this problem as I am not experienced with the beta function.
Thank you.

Note: Solution should be done without using harmonic series.
 A: All the following Mathematica commands calculate your limit, in decreasing order of time (the more naive one uses more time):
Limit[D[Gamma[a]*Gamma[b]/Gamma[a + b], {a, 2}, {b, 2}] /. {a -> 
  1/2 + x, b -> x}, x -> 0] // FunctionExpand // Expand

the above command directly calculates the limit, by choosing a path approaching $(a,b)=(1/2,0)$. It takes $32$ seconds on my machine.
D[Normal[Series[
        Gamma[a]*Gamma[b]/Gamma[a + b], {a, 1/2, 4}, {b, 0, 4}]], {a, 
       2}, {b, 2}] /. a -> 1/2 /. b -> 0 // FullSimplify // 
  Expand 

rather than calculating the limit, this one uses series expansion up to constant term. It takes $12$ seconds.
Normal[Series[
     D[Gamma[a]*Gamma[b]/Gamma[a + b], {a, 2}, {b, 2}], {a, 1/2, 
      0}, {b, 0, 0}]] // FullSimplify // Expand

this one does not even calculate derivatives, instead uses series expansion up to 4th order. It takes $3.5$ seconds.
D[Exp[Series[
        LogGamma[a] + LogGamma[b] - LogGamma[a + b], {a, 1/2, 4}, {b, 
         0, 4}]], {a, 2}, {b, 2}] /. a -> 1/2 /. b -> 0 // 
   FullSimplify // Expand

This use the well-known simple series of log gamma function, it takes only $0.5$ seconds.

It's easy to guess why the fourth one is most efficient. To see how much is used for each computation, execute ClearSystemCache[];(your command)//Timing.
Such beta limit arising from logarithm integrals is well-known, it's also not difficult to write down an recursion for it.
A: A magnificent solution by Cornel without using the derivative of beta function.

We have the identity 
$$\frac{\ln^2(1-x)}{1-x}=\sum_{n=1}^\infty (H_n^2-H_n^{(2)})x^n$$
replace $x$ with $x^2$, then multiply both sides by $\ln^2x$ and integrate from $x=0$ to $1$ we get
$$\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx=2\sum_{n=1}^\infty\frac{H_n^2-H_n^{(2)}}{(2n+1)^3}\tag1$$

By the master theorem, we have 
$$3n\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(k+1)(k+n+1)}=H_n^3+3H_nH_n^{(2)}+2H_n^{(3)}$$
differentiate both sides with respect to $n$ to get
$$3\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(k+n+1)^2}=\frac{d}{dn}\left(H_n^3+3H_nH_n^{(2)}+2H_n^{(3)}\right)$$
differentiate again and let $n\mapsto -1/2$ we get
$$-6\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(k+1/2)^3}=\frac{d^2}{dn^2}\left(H_n^3+3H_nH_n^{(2)}+2H_n^{(3)}\right)_{n\mapsto-1/2}$$
or 
$$\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(2k+1)^3}=\frac{31}{2}\zeta(5)-\frac{45}{8}\ln2\zeta(4)+\frac72\ln^22\zeta(3)-7\zeta(2)\zeta(3)\tag2$$

From (1) and (2) we get

$$\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx=31\zeta(5)-\frac{45}{4}\ln2\zeta(4)+7\ln^22\zeta(3)-14\zeta(2)\zeta(3)$$

A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
I & \equiv
{1 \over 8}\int_{0}^{1}{\ln^{2}\pars{y}\ln^{2}\pars{1 - y} \over \root{y}\pars{1 - y}}\,\dd y
\\[5mm] & =
\left.{1 \over 8}{\partial^{2} \over
\partial \mu^{2}\,\partial \nu^{2}}\int_{0}^{1}
\pars{y^{\,\mu - 1/2} - 1}\pars{1 - y}^{\,\nu - 1} \,\dd y\,
\right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}}
\\[5mm] & =
\left.{1 \over 8}{\partial^{4} \over
\partial \mu^{2}\,\partial \nu^{2}}\bracks{{\Gamma\pars{\mu +1/2}
\Gamma\pars{\nu} \over \Gamma\pars{\mu + \nu + 1/2}} -
{1 \over \nu}}\,
\right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}}
\\[5mm] & =
\left.{1 \over 8}{\partial^{4} \over
\partial \mu^{2}\,\partial \nu^{2}}\braces{{1 \over \nu}\bracks{{\Gamma\pars{\mu +1/2}
\Gamma\pars{\nu + 1} \over \Gamma\pars{\mu + \nu + 1/2}} - 1}}\,
\right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}}
\\[5mm] & =
\left.{1 \over 24}{\partial^{5} \over
\partial \mu^{2}\,\partial \nu^{3}}{\Gamma\pars{\mu +1/2}
\Gamma\pars{\nu + 1} \over \Gamma\pars{\mu + \nu + 1/2}}\,
\right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}}
\\[5mm] & =
\bbx{-\,{7 \pi ^2 \over 3}\,\zeta\pars{3} + 31\zeta\pars{5} + 7\zeta\pars{3}\ln^{2}\pars{2} - {1 \over 8}\,\pi^{4}\ln\pars{2}}
\\[5mm] & \approx 0.0654
\end{align}
