Let the considered integral be denoted by $I$. Our starting point is to reduce the number of logarithms of different arguments in the integrand. Thus, using the fact
that $6ab^2=(a+b)^3-2a^3+(a-b)^3$ we obtain
\begin{align*}
6I&=\underbrace{\int_0^1\frac{\log x}{x}\log^3(1-x^2)dx}_{x\leftarrow \sqrt{u}}
-2\int_0^1\frac{\log x}{x}\log^3(1-x)dx+
\underbrace{\int_0^1\frac{\log x}{x}\log^3\left(\frac{1-x}{1+x}\right)dx}_{x\leftarrow\frac{1-u}{1+u}}
\\
&= \frac{1}{4}\int_0^1\frac{\log u}{u}\log^3(1-u)du
-2\int_0^1\frac{\log x}{x}\log^3(1-x)dx+
2\int_0^1\frac{\log^3 u}{1-u^2}\log\left(\frac{1-u}{1+u}\right)du \\
&=\frac{-7}{4}\int_0^1\frac{\log x}{x}\log^3(1-x)dx+
2\int_0^1\frac{\log^3 x}{1-x^2}\log\left(\frac{1-x}{1+x}\right)dx \\
\end{align*}
Thus,
$$
I=-\frac{7}{24}J+\frac{1}{3}K \tag{1}
$$
with
$$
J=\int_0^1\frac{\log(1- x)}{1-x}\log^3x\,dx,\quad
K=\int_0^1\frac{1}{1-x^2}\log\left(\frac{1-x}{1+x}\right)\log^3 x\,dx
\tag{2}
$$
Now, for $x\in(-1,1)$, we have
$$
\frac{\log(1-x)}{1-x}=-\left(\sum_{k=0}^\infty x^k\right)\left(\sum_{k=1}^\infty \frac{x^k}{k}\right)
=-\sum_{n=1}^\infty H_nx^n
$$
where $H_n=\sum_{k=1}^n1/k$ is the $n$-th Harmonic number. Similarly,
$$
\frac{1}{1-x^2}\log\left(\frac{1-x}{1+x}\right)
=-\left(\sum_{k=0}^\infty x^{2k}\right)\left(\sum_{k=1}^\infty \frac{-2x^{2k-1}}{2k-1}\right)
=-2\sum_{n=1}^\infty \widetilde{H}_nx^{2n-1}
$$
where
$$\widetilde{H}_n=\sum_{k=1}^n\frac{1}{2k-1}
=H_{2n}-\frac{1}{2}H_n$$
Thus, using the fact that $\int_0^1x^n(-\log x)^3\,dx=-6/(n+1)^4$ we conclude that
\begin{align*}
J&=\sum_{n=1}^\infty H_n\int_0^1x^n(-\log x)^3\,dx=6\sum_{n=1}^\infty
\frac{H_n}{(n+1)^4}\\
&=6\sum_{n=0}^\infty
\frac{1}{(n+1)^4}\left(H_{n+1}-\frac{1}{n+1}\right)=6A-6\zeta(5)\tag{3}
\end{align*}
with
\begin{equation*}
A=\sum_{n=1}^\infty\frac{H_n}{n^4}\tag{4}
\end{equation*}
Similarly,
\begin{align*}
K&=\sum_{n=1}^\infty (2H_{2n}-H_n)\int_0^1x^{2n-1}(-\log x)^3\,dx=6\sum_{n=1}^\infty\frac{2H_{2n}-H_n}{(2n)^4}\\
&=6\sum_{n=1}^\infty\frac{(1+(-1)^{2n})H_{2n}}{(2n)^4}
-\frac{3}{8}\sum_{n=1}^\infty\frac{H_n}{n^4}\\
&=6\sum_{n=1}^\infty\frac{(1+(-1)^{n})H_{n}}{n^4}
-\frac{3}{8}\sum_{n=1}^\infty\frac{H_n}{n^4}\\
&=\frac{45}{8}\sum_{n=1}^\infty\frac{H_{n}}{n^4}
+6\sum_{n=1}^\infty\frac{(-1)^{n}H_{n}}{n^4}=\frac{45}{8}A+6B\tag{5}
\end{align*}
with
\begin{equation*}
B=\sum_{n=1}^\infty\frac{(-1)^nH_n}{n^4}\tag{6}
\end{equation*}
Combining (1) with (3) and (5) we find that
\begin{equation*}
I=\frac{1}{8}A+2B+\frac{7}{4}\zeta(5)\tag{7}
\end{equation*}
Now, sums $A$ and $B$ are known (see here ) (in a more general setting), and
we have
$$
A=3\zeta(5)-\zeta(2)\zeta(3),\qquad
B=-\frac{59}{32}\zeta(5)+\frac{1}{2}\zeta(2)\zeta(3)
$$
Thus
$$
I=\frac{7}{8}\zeta(2)\zeta(3)-\frac{25}{16}\zeta(5).
$$
which is the announced result.$\qquad\square$
I usually don't answer my own questions, but I literally just derived the solution and I think it's beautiful.
So we know from the link in the OP that:
$$\zeta(s,a)=\frac{1}{\Gamma(s)} \int_0^\infty \frac{z^{s-1} dz}{e^{a z} (1-e^{-z})}$$
Which makes our expression:
$$\zeta \left(\frac12, \frac{t}{2}\right)-\zeta \left(\frac12, \frac{t+1}{2}\right)=\frac{1}{\sqrt{\pi}} \int_0^\infty \frac{e^{- t z/2}(1-e^{-z/2}) dz}{(1-e^{-z}) \sqrt{z}}=$$
$$=\frac{2}{\sqrt{\pi}} \int_0^\infty \frac{e^{- t u^2/2}(1-e^{-u^2/2}) du}{1-e^{-u^2}}$$
Now we can see that it's very easy to take integral over $t$ (integration by parts):
$$\int_0^1 \sin \pi t~ e^{- t u^2/2} dt=\frac{4 \pi}{4 \pi^2 +u^4} (1+e^{-u^2/2})$$
Now we substitute this into the second integral to get (this is the beautiful part):
$$8 \sqrt{\pi} \int_0^\infty \frac{(1+e^{-u^2/2})(1-e^{-u^2/2}) du}{(1-e^{-u^2})(4 \pi^2 +u^4)}=8 \sqrt{\pi}\int_0^\infty \frac{du}{4 \pi^2 +u^4}$$
After a change of variables we have:
$$\frac{2 \sqrt{2}}{\pi} \int_0^\infty \frac{dv}{1 +v^4}=\frac{2 \sqrt{2}}{\pi} \frac{\pi}{2 \sqrt{2}}=1$$
Just as it was supposed to be.
God, Mathematics is simply perfect sometimes.
Appendix
Trying to prove in a simple way that the integral formula for $\zeta(1/2,a)$ equals the series definition.
$$ \int_0^\infty \frac{z^{-1/2} dz}{e^{a z} (1-e^{-z})}=2 \int_0^\infty \frac{e^{-a u^2}du}{ 1-e^{-u^2}}=2 \sum_{n=0}^\infty \int_0^\infty e^{-(a+n) u^2} du=$$
$$=2 \sum_{n=0}^\infty \frac{1}{\sqrt{a+n}} \int_0^\infty e^{-w^2} dw= \sqrt{\pi} \sum_{n=0}^\infty \frac{1}{\sqrt{a+n}} $$
Formally, this exactly fits the series definition, but it doesn't converge (it's alright, as the function for $s=1/2$ is defined by analytic continuation).
On the other hand, for the particular function in my case, the proof works well, as the series inside the integral becomes alternating due to a factor $(1-e^{-u^2/2})$ in the numerator, so everyting converges.
I would say that all of this constitutes a nice real proof of the Fresnel integrals, especially since the Poisson integral also has a few real proofs.
Though I'm definitely not claiming this is a new result. It was new for me, but a two second google search found me this paper https://www.jstor.org/stable/2320230, and I'm sure there's plenty more.
Best Answer
Here are some ideas towards explaining the form of the right hand side. I'm a bit stuck and my main approach hasn't worked out. This may just be rephrasing things in terms of other log-integrals, but hopefully this is a useful way of looking at the problem.
Taking the integral $$ I = \int_0^1 \log^2(1-x) \log^2(x) \log^3(1+x) \frac{dx}{x} $$ we can also rewrite this as $$ I = \int_0^\infty \log^2(1-e^{-x}) \log^2(e^{-x}) \log^3(1+e^{-x}) \; dx $$ which is suited for interpretation as a Mellin transform. Specifically, the power of $x$, is controlled by the power on $\log(x)$ in the original integral format as $$ I = \int_0^\infty x^2 \log^2(1-e^{-x})\log^3(1+e^{-x}) \; dx $$ according to Mathematica we have in general a result for the Mellin transform of the other components $$ \mathcal{M}[\log^n(1\pm e^{-x})](s) = (-1)^n n! \Gamma(s) S_{s,n}(\mp 1) $$ invoking the Neilsen Generalisation of the polylogarithm, $S_{s,n}$. This does recreate the series expansion for $\log(1+e^{-x})$ but the series for $\log(1-e^{-x})$ has a $\log(x)$ term, which might be causing a problem.
We could toy with the idea of a formal series via the Ramanujan Master Theorem, using these Mellin transforms $$ \log^n(1\pm e^{-x}) = \sum_{k=0}^\infty \frac{(-1)^{k+n} n!}{k!} S_{-k,n}(\mp 1)x^k $$ and then the Cauchy product $$ \log^a(1 + e^{-x})\log^b(1 - e^{-x}) = \left( \sum_{k=0}^\infty \frac{(-1)^{k+a} a!}{k!} S_{-k,a}(-1)x^k \right)\left( \sum_{k=0}^\infty \frac{(-1)^{k+b} b!}{k!} S_{-k,b}(1)x^k \right) $$ $$ \log^a(1 + e^{-x})\log^b(1 - e^{-x}) = \sum_{k=0}^\infty \left(\sum_{l=0}^k \frac{(-1)^{a+b+k} a! b!}{l!(k-l)!} S_{-l,a}(-1) S_{l-k,b}(1)\right) x^k $$ alternatively $$ \log^a(1 + e^{-x})\log^b(1 - e^{-x}) = \sum_{k=0}^\infty \frac{(-1)^k}{k!} \left(\sum_{l=0}^k (-1)^{a+b} a! b! \binom{k}{l} S_{-l,a}(-1) S_{l-k,b}(1)\right) x^k $$ plausibly leading to (via RMT) $$ \mathcal{M}\left[ \log^a(1 + e^{-x})\log^b(1 - e^{-x})\right](s) = \Gamma(s) \sum_{l=0}^{-s} (-1)^{a+b} a! b! \binom{-s}{l} S_{-l,a}(-1) S_{l-k,b}(1) $$ then we would conceptually have (with some dodgy negative parts) an answer for the integral as a sum over (four?) pairs of generalized Polylogs, specifically in the case that $s=3$.
This motivates an expression in terms of pairs of $S_{n,k}(z)$, we can guess a term and quickly find $$ -8\cdot3 \cdot 19 S_{2,2}(1)S_{1,3}(-1) = -\frac{19}{15} \pi ^4 \text{Li}_4\left(\frac{1}{2}\right)-\frac{133}{120} \pi ^4 \zeta (3) \log (2)+\frac{19 \pi ^8}{1350}+\frac{19}{360} \pi ^6 \log ^2(2)-\frac{19}{360} \pi ^4 \log ^4(2) $$ this covers a few of the terms in your expression R.H.S. It is likely that other terms contribute to $\pi^8$ for example. I can't get an explicit value for $S_{2,3}(-1)$ to explore this further, but I would assume this holds a $\mathrm{Li}_5(1/2)$ term among others, and the other factor is $S_{1,2}(1) = \zeta(3)$. Perhaps your linear combinations method can be rephrased in terms of the generalized polylogarithm?