Evaluate $\int^1_0 \log^2(1-x) \log^2(x) \, dx$

Using the series of $\ln^2(1-x)$, \begin{align} \int^1_0\ln^2{x}\ln^2(1-x) \ {\rm d}x &=\sum^\infty_{n=1}\frac{2H_n}{n+1}\int^1_0x^{n+1}\ln^2{x}\ {\rm d}x\\ &=\sum^\infty_{n=1}\frac{4H_n}{(n+1)(n+2)^3} \end{align} Then integrate $f(z)=\dfrac{(\gamma+\psi_0(-z))^2}{(z+1)(z+2)^3}$ along an infinitely large square. The integral vanishes which implies the sum of its residues is zero. At the positive integers, \begin{align} \sum^\infty_{n=1}{\rm Res}(f,n) &=\sum^\infty_{n=1}\operatorname*{Res}_{z=n}\left[\frac{1}{(z+1)(z+2)^3(z-n)^2}+\frac{2H_n}{(z+1)(z+2)^3(z-n)}\right]\\ &=\sum^\infty_{n=1}\frac{2H_n}{(n+1)(n+2)^3}-\sum^\infty_{n=1}\frac{4n+5}{(n+1)^2(n+2)^4}\\ &=\sum^\infty_{n=1}\frac{2H_n}{(n+1)(n+2)^3}+\frac{\pi^4}{30}+2\zeta(3)-\frac{91}{16} \end{align} At $z=0$, \begin{align} {\rm Res}(f,0) &=\operatorname*{Res}_{z=0}\frac{1}{z^2(z+1)(z+2)^3}\\ &=-\frac{5}{16} \end{align} At $z=-1$, \begin{align} {\rm Res}(f,-1) &=0 \end{align} At $z=-2$, \begin{align} {\rm Res}(f,-2) &=\frac{1}{2}\lim_{z\to-2}\frac{{\rm d}^2}{{\rm d}z^2}\frac{(\gamma+\psi_0(-z))^2}{z+1}\\ &=-\frac{\pi^4}{36}+2\zeta(3)+\frac{2\pi^2}{3}-6 \end{align} Therefore, \begin{align} \int^1_0\ln^2{x}\ln^2(1-x)\ {\rm d}x &=2\left[-\frac{\pi^4}{30}-2\zeta(3)+\frac{91}{16}+\frac{5}{16}+\frac{\pi^4}{36}-2\zeta(3)-\frac{2\pi^2}{3}+6\right]\\ &=-\frac{\pi^4}{90}-8\zeta(3)-\frac{4\pi^2}{3}+24 \end{align}

Starting with the Beta function \begin{align} B(x, y) = \int_{0}^{1} t^{x-1} (1-t)^{y-1} \, dt \end{align} differentiate with respect to $x$ and $y$ twice. This leads to \begin{align} I &= \int_{0}^{1} t^{x-1} (1-t)^{y-1} \, \ln^{2}(t) \, \ln^{2}(1-t) \, dt \\ &= \partial_{x}^{2} \partial_{y}^{2} B(x,y). \end{align} Now the integral in question is the case for $x=y=1$. For this, it is seen that \begin{align} \partial_{x}^{2} \left. B(x,y) \right|_{x=1} = \frac{1}{y} \left[ \psi'(1) - \psi'(y+1) + \left( \psi(1) - \psi(y+1) \right)^{2} \right]. \end{align} Differential with respect to $y$ yields \begin{align} \partial_{y}^{2} \partial_{x}^{2} \left. B(x,y) \right|_{x=1} &= \frac{1}{y} \left[ - \psi'''(y+1) - 2 \psi''(y+1) \left( \psi(1) - \psi(y+1) \right) + 2 \left( \psi'(y+1) \right)^{2} \right] \\ & \hspace{5mm} - \frac{2}{y^{2}} \left[ - \psi''(y+1) - 2 \psi'(y+1) \left( \psi(1) - \psi(y+1) \right) \right] \\ & \hspace{10mm} + \frac{2}{y^{3}} \left[ \psi'(1) - \psi'(y+1) + \left( \psi(1) - \psi(y+1) \right)^{2} \right]. \end{align} Now \begin{align} \partial_{y}^{2} \partial_{x}^{2} \left. B(x,y) \right|_{x=y=1} &= 4 - \psi'''(2) + 4 \psi''(2) - 4 \psi'(2) + 2 \left( \psi'(2) \right)^{2} . \end{align} Since \begin{align} \psi'(2) &= \frac{\pi^{2}}{6} -1 \\ \psi''(2) &= 2 - 2\zeta(3) \\ \psi'''(2) &= \frac{\pi^{4}}{45} - 6 \end{align} then \begin{align} \partial_{y}^{2} \partial_{x}^{2} \left. B(x,y) \right|_{x=y=1} = 24 - \frac{4 \pi^{2}}{3} - \frac{\pi^{4}}{90} - 8 \zeta(3). \end{align}

First, letting $x=e^{-t}$ gives $$ \begin{align} \int_0^1x^k\log(x)^2\,\mathrm{d}x &=\int_0^\infty e^{-(k+1)t}\,t^2\,\mathrm{d}t\\ &=\frac2{(k+1)^3}\tag{1} \end{align} $$ Also, squaring the power series yields $$ \log(1-x)^2=\sum_{k=1}^\infty\frac{2H_{k-1}}{k}x^k\tag{2} $$ Therefore, $$ \begin{align} &\int_0^1\log(1-x)^2\log(x)^2\,\mathrm{d}x\\ &=\int_0^1\sum_{k=1}^\infty\frac{2H_{k-1}}{k}x^k\log(x)^2\,\mathrm{d}x\tag{3a}\\ &=\sum_{k=1}^\infty\frac{4H_{k-1}}{k(k+1)^3}\tag{3b}\\ &=\sum_{k=1}^\infty\sum_{j=1}^{k-1}\frac4{jk(k+1)^3}\tag{3c}\\ &=\sum_{k=1}^\infty\sum_{j=1}^k\frac4{jk(k+1)^3}-\sum_{k=1}^\infty\frac4{k^2(k+1)^3}\tag{3d}\\ &=\sum_{j=1}^\infty\frac4j\sum_{k=j}^\infty\left(\frac1k-\frac1{k+1}-\frac1{(k+1)^2}-\frac1{(k+1)^3}\right)\tag{3e}\\ &-4\sum_{k=1}^\infty\left(-\frac3k+\frac3{k+1}+\frac1{k^2}+\frac2{(k+1)^2}+\frac1{(k+1)^3}\right)\tag{3f}\\ &=4\zeta(2)-\sum_{j=1}^\infty\frac4j\left[\sum_{k=j}^\infty\left(\frac1{k^2}+\frac1{k^3}\right)-\frac1{j^2}-\frac1{j^3}\right]\tag{3g}\\[6pt] &-4(-6+3\zeta(2)+\zeta(3))\tag{3h}\\[15pt] &=24-8\zeta(2)-4\zeta(3)\tag{3i}\\[4pt] &-\sum_{k=1}^\infty\sum_{j=1}^k\frac4j\left(\frac1{k^2}+\frac1{k^3}\right)+4\zeta(3)+4\zeta(4)\tag{3j}\\[2pt] &=24-8\zeta(2)+4\zeta(4)\tag{3k}\\[6pt] &-\sum_{k=1}^\infty\frac{4H_k}{k^2}-\sum_{k=1}^\infty\frac{4H_k}{k^3}\tag{3l}\\ &=24-8\zeta(2)-6\zeta(4)-8\zeta(3)+2\zeta(2)^2\tag{3m}\\[6pt] &=24-\frac{4\pi^2}3-\frac{\pi^4}{90}-8\zeta(3)\tag{3n} \end{align} $$ Explanation:
$\text{(3a)}$: apply $(2)$
$\text{(3b)}$: apply $(1)$
$\text{(3c)}$: expand $H_{k-1}$
$\text{(3d)}$: add and subtract the $j=k$ term
$\text{(3e)}$: partial fractions: $\frac1{k(k+1)^3}=\frac1k-\frac1{k+1}-\frac1{(k+1)^2}-\frac1{(k+1)^3}$
$\text{(3f)}$: partial fractions: $\frac1{k^2(k+1)^3}=-\frac3k+\frac3{k+1}+\frac1{k^2}+\frac2{(k+1)^2}+\frac1{(k+1)^3}$
$\text{(3g)}$: apply $\sum\limits_{j=1}^\infty\frac4j\sum\limits_{k=j}^\infty\left(\frac1k-\frac1{k+1}\right)=\sum\limits_{j=1}^\infty\frac4{j^2}=4\zeta(2)$ and
$\hphantom{\text{(3g):}}$ $\sum\limits_{k=j}^\infty\left(\frac1{(k+1)^2}+\frac1{(k+1)^3}\right)=\sum\limits_{k=j}^\infty\left(\frac1{k^2}+\frac1{k^3}\right)-\frac1{j^2}-\frac1{j^3}$ to $\text{(3e)}$
$\text{(3h)}$: evaluate $\text{(3f)}$
$\text{(3i)}$: collect the evaluated terms in $\text{(3g)}$ and $\text{(3h)}$
$\text{(3j)}$: change order of summation in $\text{(3g)}$ and evaluate sums
$\text{(3k)}$: collect the evaluated terms in $\text{(3i)}$ and $\text{(3j)}$
$\text{(3l)}$: collect $H_k$ in $\text{(3j)}$
$\text{(3m)}$: use this answer to evaluate the sums in $\text{(3l)}$
$\text{(3n)}$: use $\zeta(2)=\frac{\pi^2}6$ and $\zeta(4)=\frac{\pi^4}{90}$

(Too long for comment)

For what it's worth, I'd like to expand on my comment above about attacking this integral via integration by parts. My main reason for wanting to solve the problem this way is not because it is the easiest way and certainly not because it is the quickest way, but rather because it seems to be the most elementary way, relying for the most part on methods accessible to high school students. (My apologies to the OP for just now getting around to writing this, but partially due to the length of the derivation I simply didn't have the time until now.)

We may reduce the evaluation of the definite integral $\mathcal{I}:=\int_{0}^{1}\ln^2{(x)}\ln^2{(1-x)}\,\mathrm{d}x$ to the problem of evaluating a sum of relatively less difficult integrals via integration by parts (IBPs). We'll use $u$ and $u^\prime$ to denote the following functions:

$$\begin{cases} u{(x)}:=\ln^2{(x)}\ln{(1-x)},\\ v^{\prime}{(x)}:=\ln{(1-x)}.\\ \end{cases}$$

The derivative of $\ln^2{(x)}\ln{(1-x)}$ is:

$$\frac{d}{dx}\left(\ln^2{(x)}\ln{(1-x)}\right)=\frac{2\ln{(x)}\ln{(1-x)}}{x}-\frac{\ln^2{(x)}}{1-x};$$

the indefinite integral of $\ln{(1-x)}$ is:

$$\begin{align} \int\ln{(1-x)}\,\mathrm{d}x &=x\ln{(1-x)}-\int\frac{x}{x-1}\,\mathrm{d}x\\ &=x\ln{(1-x)}-x-\ln{(1-x)}+\color{grey}{constant}. \end{align}$$

Then, choosing $v{(x)}$ by selecting a null constant term in the last line above, we have:

$$\begin{cases} u^{\prime}{(x)}=\frac{2\ln{(x)}\ln{(1-x)}}{x}-\frac{\ln^2{(x)}}{1-x},\\ v{(x)}=x\ln{(1-x)}-x-\ln{(1-x)}.\\ \end{cases}$$

IBPs in this way yields:

$$\begin{align} \mathcal{I} &=\int_{0}^{1}\ln^2{(x)}\ln^2{(1-x)}\,\mathrm{d}x\\ &=\int_{0}^{1}u{(x)}\,v^{\prime}{(x)}\,\mathrm{d}x\\ &=\left[u{(x)}\,v{(x)}\right]_{0}^{1}-\int_{0}^{1}u^{\prime}{(x)}\,v{(x)}\,\mathrm{d}x\\ &=\left[\ln^2{(x)}\ln{(1-x)}\left(x\ln{(1-x)}-x-\ln{(1-x)}\right)\right]_{0}^{1}\\ &~~~~~ -\int_{0}^{1}\left(\frac{2\ln{(x)}\ln{(1-x)}}{x}-\frac{\ln^2{(x)}}{1-x}\right)\left(x\ln{(1-x)}-x-\ln{(1-x)}\right)\,\mathrm{d}x.\\ \end{align}$$

Note that after integrating by parts both boundary terms vanish identically on the interval $[0,1]$. Thus,

$$\begin{align} \mathcal{I} &=\left[\ln^2{(x)}\ln{(1-x)}\left(x\ln{(1-x)}-x-\ln{(1-x)}\right)\right]_{0}^{1}\\ &~~~~~ -\int_{0}^{1}\left(\frac{2\ln{(x)}\ln{(1-x)}}{x}-\frac{\ln^2{(x)}}{1-x}\right)\left(x\ln{(1-x)}-x-\ln{(1-x)}\right)\,\mathrm{d}x\\ &=0-0+\int_{0}^{1}\left(\frac{\ln^2{(x)}}{1-x}-\frac{2\ln{(x)}\ln{(1-x)}}{x}\right)\left(x\ln{(1-x)}-x-\ln{(1-x)}\right)\,\mathrm{d}x\\ &=\int_{0}^{1}\left(\frac{\ln^2{(x)}}{1-x}-\frac{2\ln{(x)}\ln{(1-x)}}{x}\right)\left(x\ln{(1-x)}-x-\ln{(1-x)}\right)\,\mathrm{d}x\\ &=\color{green}{\int_{0}^{1}\frac{x\ln^2{(x)}\ln{(1-x)}}{1-x}\,\mathrm{d}x}-\int_{0}^{1}\frac{x\ln^2{(x)}}{1-x}\,\mathrm{d}x\color{green}{-\int_{0}^{1}\frac{\ln^2{(x)}\ln{(1-x)}}{1-x}\,\mathrm{d}x}\\ &~~~~~ \color{brown}{-\int_{0}^{1}2\ln{(x)}\ln^2{(1-x)}\,\mathrm{d}x}+\int_{0}^{1}2\ln{(x)}\ln{(1-x)}\,\mathrm{d}x\color{brown}{+\int_{0}^{1}\frac{2\ln{(x)}\ln^2{(1-x)}}{x}\,\mathrm{d}x}\\ &=\color{green}{-\int_{0}^{1}\ln^2{(x)}\ln{(1-x)}\,\mathrm{d}x}\color{orange}{-\int_{0}^{1}\frac{x\ln^2{(x)}}{1-x}\,\mathrm{d}x}+2\int_{0}^{1}\ln{(x)}\ln{(1-x)}\,\mathrm{d}x\\ &~~~~~ \color{brown}{+2\int_{0}^{1}\frac{(1-x)\ln{(x)}\ln^2{(1-x)}}{x}\,\mathrm{d}x}\\ &=2\int_{0}^{1}\ln{(x)}\ln{(1-x)}\,\mathrm{d}x\color{orange}{+\int_{0}^{1}\mathrm{d}x-\int_{0}^{1}\frac{\ln^2{(x)}}{1-x}\,\mathrm{d}x}-\int_{0}^{1}\ln^2{(x)}\ln{(1-x)}\,\mathrm{d}x\\ &~~~~~ +2\int_{0}^{1}\frac{(1-x)\ln{(x)}\ln^2{(1-x)}}{x}\,\mathrm{d}x\\ &=1+2\int_{0}^{1}\ln{(x)}\ln{(1-x)}\,\mathrm{d}x-\int_{0}^{1}\frac{\ln^2{(x)}}{1-x}\,\mathrm{d}x-\int_{0}^{1}\ln^2{(x)}\ln{(1-x)}\,\mathrm{d}x\\ &~~~~~ +2\int_{0}^{1}\frac{(1-x)\ln{(x)}\ln^2{(1-x)}}{x}\,\mathrm{d}x\\ &=:1+2\,\mathcal{I}_{1}-\mathcal{I}_{2}-\mathcal{I}_{3}+2\,\mathcal{I}_{4}, \end{align}$$

where in the last line we've simply introduced the constant names $\mathcal{I}_{k},k\in\{1,2,3,4\}$, for convenience to denote the following integrals, respectively:

$$\begin{cases} \int_{0}^{1}\ln{(x)}\ln{(1-x)}\,\mathrm{d}x,\\ \int_{0}^{1}\frac{\ln^2{(x)}}{1-x}\,\mathrm{d}x,\\ \int_{0}^{1}\ln^2{(x)}\ln{(1-x)}\,\mathrm{d}x,\\ \int_{0}^{1}\frac{(1-x)\ln{(x)}\ln^2{(1-x)}}{x}\,\mathrm{d}x.\\ \end{cases}$$

Now, the OP noted in the question statement that WolframAlpha was unable to find a closed form for the definite integral $\mathcal{I}$. As one might expect, WA is also unable to find an anti-derivative for the corresponding indefinite integral. HOWEVER, it turns out that WolframAlpha does return anti-derivatives for each of the four integrals $\mathcal{I}_{k},k\in\{1,2,3,4\}$. WA is also capable of returning closed form values for each of the definite integrals except $\mathcal{I}_{4}$, which has the longest and most complicated anti-derivative of the bunch.