A second solution in large steps by Cornel Ioan Valean
Let's start with the following useful identity which is easily derived by using recurrence relations and simple rearrangements, manipulations with sums, that is
Let $n$ be a non-negative integer number. Then, we have
$$\int_0^1 x^{2n}\frac{\log(1+x)}{1+x}\textrm{d}x$$
$$=\frac{1}{2}H_{2n}^2-2\log(2) H_{2n}+\frac{1}{2}H_{2n}^{(2)}-\frac{1}{4}H_n^2-\frac{1}{4}H_n^{(2)}+\log (2)H_n+\frac{1}{2} \log ^2(2)-\sum_{k=1}^{n-1}\frac{H_k}{2 k+1},$$
where $H_n^{(m)}=1+\frac{1}{2^m}+\cdots+\frac{1}{n^m}$ represents the $n$th generalized harmonic number of order $m$.
By multiplying both sides of the identity above by $1/n^3$ and considering the summation from $n=1$ to $\infty$, we get
$$\sum_{n=1}^{\infty} \frac{1}{n^3}\sum_{k=1}^{n-1}\frac{H_{k}}{2 k+1}=\sum_{k=1}^{\infty} \sum_{n=k+1}^{\infty}\frac{1}{n^3}\frac{H_{k}}{2 k+1}=\underbrace{\sum_{k=1}^{\infty}\frac{H_{k}}{2 k+1}\left(\zeta(3)-H_k^{(3)}\right)}_{\text{The desired series}}$$
$$=\frac{1}{2}\sum_{n=1}^{\infty}\frac{H_{2n}^2}{n^3}-2\log(2) \sum_{n=1}^{\infty}\frac{H_{2n}}{n^3}+\frac{1}{2}\sum_{n=1}^{\infty}\frac{H_{2n}^{(2)}}{n^3}-\frac{1}{4}\sum_{n=1}^{\infty}\frac{H_n^2}{n^3}-\frac{1}{4}\sum_{n=1}^{\infty} \frac{H_n^{(2)}}{n^3}$$
$$+\log (2)\sum_{n=1}^{\infty} \frac{H_n}{n^3}+\frac{1}{2}\log ^2(2)\sum_{n=1}^{\infty}\frac{1}{n^3}-\int_0^1 \frac{\log(1+x)}{1+x}\operatorname{Li}_3(x^2)\textrm{d}x,$$
where we see all the series in the right-hand side are easily reducible to known series which may also be found in the book (Almost) Impossible Integrals, Sums, and Series.
On the other hand, with simple integration by parts, we obtain
$$\int_0^1 \frac{\log(1+x)}{1+x}\operatorname{Li}_3(x^2)\textrm{d}x$$
$$=\frac{1}{2}\log^2(2)\zeta(3)-2\int_0^1 \frac{\log^2(1+x)\operatorname{Li}_2(x)}{x}\textrm{d}x-2\int_0^1 \frac{\log^2(1+x)\operatorname{Li}_2(-x)}{x}\textrm{d}x,$$
where the last integrals may be found calculated in the paper The calculation of a harmonic series with a weight $5$ structure, involving the product of harmonic numbers, $H_n H_{2n}^{(2)}$.
A note: The sister of the result above (easy to obtain by recurrence relations and very useful),
$$\int_0^1 x^{2n-1} \frac{\log(1+x)}{1+x}\textrm{d}x$$
$$=2\log(2) H_{2n}-\log(2)H_n+\frac{1}{4}H_n^2+\frac{1}{4}H_n^{(2)}-\frac{1}{2}H_{2n}^2-\frac{1}{2} H_{2n}^{(2)}+\frac{H_{2n}}{2n}-\frac{H_n}{2n}
$$
$$
-\frac{1}{2}\log^2(2)+\sum_{k=1}^{n-1}\frac{H_k}{2 k+1}.
$$
Using the Pochhammer symbol notation, let $$S_0 = \sum_{n=0}^\infty \frac{(1/2)_n^3}{n!^3}\sum_{k=1}^n \frac{1}{2k-1} \qquad S_1 = \sum_{n=0}^\infty \frac{(1/2)_n^3}{n!^3}H_n$$
with $H_n$ the Harmonic number. OP's sum is $S_0$. I will show
$$\tag{0}S_0 = \frac{\Gamma \left(\frac{1}{4}\right)^4}{24 \pi ^2} \qquad S_1=\frac{(\pi -3 \log (2)) \Gamma \left(\frac{1}{4}\right)^4}{6 \pi ^3}$$
(I realized after posting this answer, that all details presented below were in fact very succinctly presented here).
First note that derivative of Pochhammer symbol can be expressed as digamma function: ${\left. {\frac{d}{{d\varepsilon }}} \right|_{\varepsilon = 0}}{(a + \varepsilon )_n} = {(a)_n}\left[ {\psi (a + n) - \psi (a)} \right]$. Consider $$\begin{aligned}&\quad {\left. {\frac{d}{{d\varepsilon }}} \right|_{\varepsilon = 0}}{_3F_2}\left( {\frac{1}{2} + \varepsilon ,\frac{1}{2} + \varepsilon ,\frac{1}{2};1 + \varepsilon ,1;1} \right)\\
&= \sum\limits_{n = 0}^\infty {\frac{{{(1/2)}_n^3}}{{n{!^3}}}\left[ {2\psi (n + \frac{1}{2}) - 2\psi (\frac{1}{2}) - \psi (1 + n) + \psi (1)} \right]} = 4S_0 - S_1
\end{aligned}$$
On other other hand, the above $_3F_2$ can be evaluated to be $$\frac{2^{-\varepsilon-\frac{1}{2}} \Gamma \left(\frac{1}{4}-\frac{\varepsilon}{2}\right) \Gamma (\varepsilon+1)}{\Gamma \left(\frac{3}{4}-\frac{\varepsilon}{2}\right) \Gamma \left(\frac{\varepsilon}{2}+\frac{3}{4}\right)^2}$$
calculating its derivative gives $$\tag{1}4S_0 - S_1 = \frac{\log (2) \Gamma \left(\frac{1}{4}\right)^4}{2 \pi ^3}$$
We need to find a second equation relating $S_0,S_1$, this is more difficult. I will state the following result, it will be proved at the end.
For $0<a<1,0<x\leq 1/2$, we have
$$\tag{*}{_2F_1}\left( {a,1 - a;1;x} \right){_2F_1}\left( {a,1 - a;1;1 - x} \right) = \frac{{\sin a\pi }}{\pi }\sum\limits_{n = 0}^\infty {{c_n}\left[ {{d_n} - \log (4x(1 - x))} \right]{{(4x(1 - x))}^n}} $$
where $${c_n} = \frac{{{{(1/2)}_n}{{(a)}_n}{{(1 - a)}_n}}}{{n{!^3}}}\qquad {d_n} = 3\psi (n + 1) - \psi (n + \frac{1}{2}) - \psi (n + a) - \psi (n + 1 - a)$$
Let $a=x=1/2$,
$$\begin{aligned}\pi {_2F_1}{\left( {\frac{1}{2},\frac{1}{2};1;\frac{1}{2}} \right)^2} &= \sum\limits_{n = 0}^\infty {\frac{{{(1/2)}_n^3}}{{n{!^3}}}\left[ {3\psi (n + 1) - 3\psi (n + \frac{1}{2})} \right]} \\
&= 3S_1 - 6S_1 +6\log 2 \sum\limits_{n = 0}^\infty {\frac{{{(1/2)}_n^3}}{{n{!^3}}}}
\end{aligned}$$
the $_3F_2$ and $_2F_1$ appear at both sides are easy (for example, an elliptic singular value for LHS, Dixon's theorem for RHS), so
$$\tag{2}3S_1 - 6S_0 = \frac{\Gamma \left(\frac{1}{4}\right)^4}{4 \pi ^2}-\frac{\Gamma \left(\frac{1}{4}\right)^4}{4 \pi ^3}(6\log 2)$$
Solving $(1), (2)$ gives $(0)$.
I briefly mention a proof of (*). For $a=1/2$, this was discovered by Watson in 1908, the following proof works essentially uniformly for all $0<a<1$.
An one-line proof: Verify both sides of $(*)$ satisfy the following 3rd order ODE:
$$(-4 a^2 x^2+4 a^2 x+4 a x^2-4 a x+6 x^2-6 x+1) y'(x)-2 (a-1) a (2 x-1) y(x)+3 (x-1) (2 x-1) x y''(x)+(x-1)^2 x^2 y^{(3)}(x)=0$$
and their first few terms of series expansion at $x=0$ agree.
A conceptual proof: $y_1(x) = {_2F_1}(a,1-a;1;x)$ satisfies the following ODE: $$(1-a) a y(x)+(2 x-1) y'(x)+(x-1) x y''(x)=0$$ which is invariant under $x\mapsto 1-x$, so $y_2(x) = y_1(1-x)$ is the other solution. This ODE is Fuchsian with a logarithmic singularity at $x=0$. On the other hand, $y_1^2 = {_3F_2}(a,1-a,1/2;1,1;4x(1-x))$. Denote $X=4x(1-x)$. $y_1^2, y_1y_2, y_2^2$ satisfies a Fuchsian 3rd order ODE in $X$. Theory of Fuchsian equation says, if $y_1(x)^2 = \sum c_n X^n$, then $$y_1(x)y_2(x) = A \log X\sum c_n X^n + B \sum d_nX^n$$
$$y_2(x)^2 = C \log^2 X \sum c_n X^n + D \log X \sum d_n X^n + E \sum e_n X^n$$
here $d_n$ (resp. $e_n$) can be explicitly given, via Frobenius method with integral exponent differences, to be $c_n$ times digamma (resp. trigamma) factors. This explains the form of $(*)$, making all these explicit is not difficult.
Best Answer
You don’t need any integrals.
$$S=\sum_i \sum_{j<=i} 1/(i^2 j^2)= \sum_i \sum_j 1/(2i^2 j^2) + \sum_i 1/(2i^4) =(\zeta(2)^2 +\zeta(4))/2= (\pi^4/36+\pi^4/90)/2=7\pi^4/360$$