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}$$