Evaluation of double summation

Here an alternative way:

$$\sum_{m\geq 2}\sum_{n=1}^{m-1}\frac{1}{(2m-1)(2n-1)(2m-2n)}=\frac{1}{2}\sum_{T\geq 1}\sum_{S=0}^{T-1}\frac{1}{(2T+1)(2S+1)(T-S)} $$ equals $$ \frac{1}{2}\sum_{S\geq 0}\sum_{T>S}\frac{1}{(2T+1)(2S+1)(T-S)}=\frac{1}{2}\sum_{S\geq 0}\sum_{d\geq 1}\frac{1}{(2S+2d+1)(2S+1)d} $$ or $$ \frac{1}{2}\sum_{S\geq 0}\sum_{d\geq 1}\int_{0}^{1}\frac{x^{2S}}{2S+1}\cdot\frac{x^{2d}}{d}\,dx=\frac{1}{2}\int_{0}^{1}\frac{\text{arctanh}(x)}{x}\cdot\left(-\log(1-x^2)\right)\,dx $$ or $$ \frac{1}{4}\int_{0}^{1}\frac{\log^2(1+x)-\log^2(1-x)}{x}\,dx $$ where $$ \int_{0}^{1}\frac{\log^2(1-x)}{x}\,dx = \int_{0}^{1}\frac{\log^2(x)}{1-x}\,dx = 2\zeta(3) $$ and $$ \int_{0}^{1}\frac{\log^2(1+x)}{x}\,dx\stackrel{\text{IBP}}{=}-2\int_{0}^{1}\frac{\log(1+x)\log(x)}{1+x}\,dx $$ is easily solved as $\frac{1}{4}\zeta(3)$ via $\log(1+x)=\log(1-x)-2\text{arctanh}(x)$, $x\mapsto\frac{1-x}{1+x}$ and

$$ \int_{0}^{1}\frac{\text{arctanh}(x)^k}{x}\,dx =\int_{0}^{+\infty}\frac{2t^k}{\sinh(2t)}\,dt=\frac{2k!}{2^k}\sum_{n\geq 0}\frac{1}{(2n+1)^{k+1}}=\frac{2k!}{2^k}\zeta(k+1)\left(1-\frac{1}{2^{k+1}}\right).$$

$$\small\sum_{m\geq 1}\frac{H_{2m-1}}{(2m+1)^2}=\sum_{m\geq 1}\frac{H_{2m}}{(2m+1)^2}-\sum_{m=1}^{+\infty}\frac{1}{2m(2m+1)^2}\stackrel{\text{PFD}}{=}\sum_{m\geq 1}\frac{H_{2m}}{(2m+1)^2}-\left(2-\frac{\pi^2}{8}-\log 2\right) $$ and in a similar way $$ \sum_{m\geq 1}\frac{H_{2m-1}}{(2m+1)^2}=\sum_{m\geq 1}\frac{H_{2m+1}}{(2m+1)^2}-\left(1-\frac{\pi^2}{8}-\log 2+\frac{7}{8}\zeta(3)\right). $$ The idea now is to evaluate $\sum_{m\geq 1}\frac{H_{2m+1}}{(2m+1)^2}$ from $\sum_{n\geq 1}\frac{H_n}{n^2}z^n$. Since $$ \sum_{n\geq 1}\frac{H_n}{n}z^n = \frac{1}{2}\log^2(1-x)+\text{Li}_2(x) $$ we have $$ \sum_{n\geq 1}\frac{H_n}{n^2}=\int_{0}^{1}\left(\frac{1}{2}\log^2(1-x)+\text{Li}_2(x)\right)\frac{dx}{x}=2\zeta(3) $$ $$ \sum_{n\geq 1}(-1)^n\frac{H_n}{n^2}=\int_{0}^{1}\left(\frac{1}{2}\log^2(1+x)+\text{Li}_2(-x)\right)\frac{dx}{x}=-\frac{5}{8}\zeta(3) $$ so $$ \sum_{k\geq 0}\frac{H_{2k+1}}{(2k+1)^2}=\frac{21}{16}\zeta(3) $$ and $$ \sum_{m\geq 1}\frac{H_{2m-1}}{(2m+1)^2} = \frac{7}{16}\zeta(3)-2+\frac{\pi^2}{8}+\log 2.$$ Anyway

$$\small \sum_{m=2}^{+\infty}\sum_{n=1}^{m-1}\frac{1}{(2m-1)(2n-1)(2m-2n)}=\sum_{m\geq 2}\frac{H_{2m-2}}{(2m-1)^2}=\sum_{m\geq 1}\frac{H_{2m}}{(2m+1)^2}=\color{red}{\frac{7}{16}\zeta(3)}. $$

Here is a shorter way. The double sum

$$s=\frac{1}{2} \sum _{m=2}^{\infty } \sum _{n=1}^{m-1} \frac{1}{(2 m-1) (2 n-1) (m-n)}\tag{1}$$

is easily done directly in integral form by using the generating function, writing

$$s = \int_{[0,1]^3}g(x,y,z) \,dx\,dy\,dz\tag{2}$$

where

$$g(x,y,z)=\frac{1}{2} \sum _{m=2}^{\infty } \sum _{n=1}^{m-1} y^{2 m-2} z^{2 n-2} x^{m-n-1}\\ =\frac{y^2}{2 \left(x y^2-1\right) \left(y^2 z^2-1\right)}\tag{3}$$

and then doing the triple integral step by step

$$g_x = \int_0^1 g \,dx = \frac{\log \left(1-y^2\right)}{2 y^2 z^2-2}\tag{4a}$$

$$g_{xz} = \int_0^1 g_x \,dz=-\frac{\log \left(1-y^2\right) \tanh ^{-1}(y)}{2 y}\\= \frac{\log ^2(1-y)-\log ^2(y+1)}{4 y} \tag{4b}$$

and finally

$$s=\int_0^1 \frac{\log ^2(1-y)-\log ^2(y+1)}{4 y} \, dy= \frac{7 \zeta (3)}{16}\tag{4c}$$

Discussion

§1. The order of integration was a lucky choice. If I would have tried the $y$-integration second that would have resulted in a complicated expression with derivatives of hypergeometric functions, and I surely would have given up instead if trying the remaining $z$-integration.

§2. A variation

For the sum with $(n-m)$ replaced by $(n+m)$ the same method takes us on a long journey through complicated structures but in the end the result is surprisingly simple:

$$s_{+}=\frac{1}{2} \sum _{m=2}^{\infty } \sum _{n=1}^{m-1} \frac{1}{(2 m-1) (2 n-1) (m+n)}\\ =\frac{3 \log (2)}{4}-\frac{\pi ^2}{32} \simeq 0.211435\tag{5}$$