Let's consider the double sum over all $\,(n,k)\in\mathbb{Z}^2\,$ except the origin $(0,0)$ :
$$\tag{1}S(s):=\sum_{(n,k)\neq (0,0)}\frac{1}{\left ( n^2+k^2 \right )^s}=4\left(\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^s}+\zeta(2\,s)\right)$$
This is a "Lorenz–Hardy sum" (ref. Lorenz L. $1871$ "Bidrag tiltalienes theori" and Hardy G. H. $1919$ "On some definite integral considered by Mellin" in his Collected papers vol $7$).
Analytical method
The Mellin transform of a function $f$ is defined by :
$$\tag{2}\{\mathcal{M}f\}(s):=\int_0^\infty t^{s-1}f(t)\;dt$$
applied to $\;f:t\mapsto e^{\large{-mt}}\;$ and from the definition of the Gamma function gives :
$$\tag{3}\frac{\Gamma(s)}{\left(m \right)^s}=\int_0^\infty t^{s-1}\,e^{-m\,t}\;dt$$
Supposing that $\Re(s)>1$ we may rewrite our double-sum as :
\begin{align}
\Gamma(s)\,S(s)&=\sum_{(n,k)\neq (0,0)}\frac{\Gamma(s)}{\left ( n^2+k^2 \right )^s}\\
&=\sum_{(n,k)\neq (0,0)}\int_0^\infty t^{s-1}\,e^{-(n^2+k^2)\,t}\;dt\\
&=\int_0^\infty t^{s-1}\sum_{(n,k)\neq (0,0)}\,e^{-(n^2+k^2)\,t}\;dt\\
&=\int_0^\infty t^{s-1}\left(\sum_{n\in\mathbb{Z}}\,e^{-n^2\,t}\sum_{k\in\mathbb{Z}}\,e^{-k^2\,t}-1\right)\;dt\\
\tag{4}&=\int_0^\infty t^{s-1}\left(\theta_3(0,e^{-t})^2-1\right)\;dt\\
\end{align}
i.e. the Mellin transform of $\,f:t\mapsto \theta_3(0,e^{-t})^2-1\,$ where the Jacobi theta function $\theta_3$ is defined by (in this answer we will implicitly suppose that $z=0$ and $q=e^{-t}$) :
$$\tag{5}\theta_3(z,q):=\sum_{n=-\infty}^\infty q^{n^2}e^{2niz}$$
In $1829$ Jacobi published his deep elliptic functions book "Fundamenta nova theoriae functionum ellipticarum" where he obtained numerous identities including the formula $4$) of page $103$ :
$$\tag{6}\theta_3(0,q)^2=\frac {2\,K}{\pi}=1+4\sum_{n=1}^\infty \frac{q^n}{1+q^{2n}}$$
($K$ is the Complete elliptic integral of the first kind but we won't use it here)
The series at the right of $(6)$ is a Lambert series. The expansion of the denominator in power series and substitution of the double-sum in $(4)$ gives :
\begin{align}
\Gamma(s)\,S(s)&=4\int_0^\infty t^{s-1}\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m q^n\,q^{2mn}\;dt\\
&=4\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m \int_0^\infty t^{s-1}e^{-n(2m+1)t}\;dt\\
&=4\,\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m \frac{\Gamma(s)}{(n(2m+1))^s}\\
&=4\,\Gamma(s)\sum_{n=1}^\infty \frac 1{n^s}\sum_{m=0}^\infty \frac{(-1)^m}{(2m+1)^s}\\
&=4\,\Gamma(s)\,\zeta(s)\,\beta(s)\\
\\
&\text{so that }\\
\\
\tag{7}\sum_{(n,k)\neq (0,0)}\frac{1}{\left ( n^2+k^2 \right )^s}&=4\,\zeta(s)\,\beta(s),\quad\Re(s)>1\\
&\text{and}\\
\tag{8}\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^s}&=\zeta(s)\,\beta(s)-\zeta(2s),\quad\Re(s)>1\\
\\\end{align}
$\qquad\qquad\qquad$with $\beta$ the Dirichlet beta function.
This powerful method allows to deduce many closed forms of lattice sums from the known theta identities.
I'll try to add one (or more) alternative derivations when time allows...
References :
- Glasser and Zucker 1980 "Lattice sums" In "Theoretical Chemistry, Advances and Perspectives Vol $5$".
- Borwein and Borwein 1987 "Pi and the AGM".
- Borwein, Glasser, McPhedran, Wan and Zucker 2013 "Lattice Sums Then and Now" (and older papers from these authors)
We obtain
\begin{align*}
\sum_{k=0}^\infty\frac{k!}{(2k+3)!!}&=\sum_{k=0}^\infty\frac{k!(2k+2)!!}{(2k+3)!}\\
&=\sum_{k=0}^\infty\frac{k!2^{k+1}(k+1)!}{(2k+3)!}\\
&=\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^k}{(2k+1)(2k+3)}\\
&=\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-1}}{2k+1}-\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-1}}{2k+3}\tag{1}
\end{align*}
We use a representation of reciprocal binomial coefficients via the Beta function:
\begin{align*}
\binom{n}{k}^{-1}=(n+1)\int_0^1z^k(1-z)^{n-k}\,dz\tag{2}
\end{align*}
and the left-hand series of (1) can be calculated as
\begin{align*}
\color{blue}{\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-1}}{2k+1}}
&=\sum_{k=0}^\infty 2^{k-1}\int_0^1z^k(1-z)^k\,dz\tag{3}\\
&=\frac{1}{2}\int_{0}^{1}\sum_{k=0}^\infty \left(2z(1-z)\right)^k\,dz\\
&=\frac{1}{2}\int_{0}^1\frac{dz}{1-2z(1-z)}\tag{4}\\
&=\frac{1}{2}\int_{0}^{1}\frac{dz}{z^2+(1-z)^2}\\
&=\frac{1}{2}\int_{0}^{\infty}\frac{du}{1+u^2}\tag{5}\\
&\,\,\color{blue}{=\frac{\pi}{4}}\tag{6}
\end{align*}
Comment:
In (3) we use the identity (2).
In (4) we apply the geometric series expansion.
In (5) we use the substitution $u=\frac{1-z}{z}, du=-\frac{1}{z^2}dz$.
We also want to apply (2) to the right-hand series of (1). To do this conveniently we need some preparatory work:
\begin{align*}
\sum_{k=0}^\infty&\binom{2k}{k}^{-1}\frac{2^{k-1}}{2k+3}\\
&=\sum_{k=0}^\infty\frac{k!k!}{(2k)!}\cdot\frac{2^{k-1}}{2k+3}\\
&=\sum_{k=0}^\infty\frac{k!(k+1)!(2k+1)}{(2k+1)!(k+1)}\cdot\frac{2^{k-1}}{2k+3}\\
&=\sum_{k=0}^\infty\frac{(k+1)!(k+1)!}{(2k+2)!}\cdot\frac{2^{k+1}}{2k+3}-\sum_{k=0}^\infty\frac{k!(k+1)!}{(2k+1)!(k+1)}\cdot\frac{2^{k-1}}{2k+3}\\
&=\sum_{k=0}^\infty\binom{2k+2}{k+1}^{-1}\frac{2^{k+1}}{2k+3}-\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-1}}{(2k+1)(2k+3)}\\
&=\sum_{k=0}^\infty\binom{2k+2}{k+1}^{-1}\frac{2^{k+1}}{2k+3}
-\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-2}}{2k+1}
+\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-2}}{2k+3}\tag{7}\\
\end{align*}
In the last line (7) we use a partial fraction decomposition as we did in (1).
We are now well prepared to do the calculation. We obtain together with (6):
\begin{align*}
\color{blue}{\sum_{k=0}^\infty \binom{2k}{k}^{-1}\frac{2^{k-2}}{2k+3}}
&=\sum_{k=0}^\infty\binom{2k+2}{k+1}^{-1}\frac{2^{k+1}}{2k+3}-\frac{\pi}{8}\\
&=\sum_{k=0}^\infty2^{k+1}\int_{0}^1z^{k+1}(1-z)^{k+1}\,dz-\frac{\pi}{8}\\
&=\sum_{k=1}^\infty2^k\int_{0}^1z^k(1-z)^k\,dz-\frac{\pi}{8}\\
&=\frac{\pi}{2}-2^0\int_{0}^1\,dz-\frac{\pi}{8}\\
&\,\,\color{blue}{=\frac{3}{8}\pi-1}\tag{8}
\end{align*}
We finally conclude from (1) together with (6) and (8)
\begin{align*}
\color{blue}{\sum_{k=0}^\infty\frac{k!}{(2k+3)!!}}
&=\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-1}}{2k+1}-\sum_{k=0}^\infty\binom{2k}{k}^{-1}\frac{2^{k-1}}{2k+3}\\
&=\frac{\pi}{4}-2\left(\frac{3}{8}\pi-1\right)\\
&\,\,\color{blue}{=2-\frac{\pi}{2}}
\end{align*}
and the claim follows.
Best Answer
Note that $$\frac{m^2+n^2}{m n\left(m^2-n^2\right)^2}= \frac{1}{2 mn(m+n)^2}+ \frac{1}{2 mn(m-n)^2}.$$ Now consider the Tornheim double sums: $$T(a,b,c)= \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} \frac{1}{m^an^b(m+n)^c}.$$ Then $$ \sum_{m,n=1\, m\neq n}^{\infty} \frac{1}{mn(m+n)^2}=T(1,1,2)-\sum_{m=1}^{\infty}\frac{1}{4m^4}=T(1,1,2)-\frac{\zeta(4)}{4},$$ $$\begin{align} \sum_{m=1}^{\infty} \sum_{n=1}^{m-1} \frac{1}{mn(m-n)^2} &= \sum_{n=1}^{\infty} \sum_{m=n+1}^{\infty} \frac{1}{mn(m-n)^2}\\ &= \sum_{n=1}^{\infty} \sum_{k=1}^{\infty} \frac{1}{nk^2(n+k)}=T(1,2,1)\end{align}$$ and again $$ \sum_{m=1}^{\infty} \sum_{n=m+1}^{\infty} \frac{1}{mn(m-n)^2}=T(1,2,1).$$ Hence $$\begin{align}\sum_{m,n=1\, m\neq n}^{\infty}{\frac{m^2+n^2}{mn(m^2-n^2)^2}} &=\frac{T(1,1,2)-\zeta(4)/4}{2} +T(1,2,1)\\ &=\frac{11\zeta(4)}{8} \end{align}$$ where we used $$T(1,1,2)=\zeta(4)/2\quad,\quad T(1,2,1)=5\zeta(4)/2$$ (see page 31 in The evaluation of Tornheim double sums).