Double summation involving factorials

Note that by letting $k=m+n$, we have that (terms can be arranged because they are non-negative), $$\begin{align} \sum^{\infty}_{m=1}\sum^{\infty}_{n=1}\frac{m\cdot n}{(m+n)!} &=\sum^{\infty}_{k=2}\frac{1}{k!}\sum_{n=1}^{k-1}n(k-n)= \frac{1}{6}\sum^{\infty}_{k=2}\frac{k^3-k}{k!}\\&=\frac{1}{6}\sum^{\infty}_{k=2}\frac{k(k-1)(k-2)+3k(k-1)}{k!}\\ &= \frac{1}{6}\sum^{\infty}_{k=3}\frac{1}{(k-3)!}+\frac{1}{2}\sum^{\infty}_{k=2}\frac{1}{(k-2)!}\\ &=\frac{e}{6}+\frac{e}{2}=\frac{2e}{3}.\end{align}$$

The slickest approaches have already been outlined by Marco and Robert, so I will show an overkill just for fun. For any $n\in\mathbb{N}$ we have $\frac{1}{n!}=[z^n]e^z = \frac{1}{2\pi i}\oint \frac{e^z}{z^{n+1}}\,dz $, hence

$$ \sum_{m,n\geq 1}\frac{mn}{(m+n)!}=\frac{1}{2\pi i}\oint\sum_{m,n\geq 1}\frac{m}{z^m}\cdot\frac{n}{z^n}\cdot\frac{e^z}{z}\,dz=\frac{1}{2\pi i}\oint\frac{z e^{z}}{(1-z)^4}\,dz $$ and $$ \sum_{m,n\geq 1}\frac{mn}{(m+n)!}=\operatorname*{Res}_{z=1}\frac{z e^z}{(1-z)^4} =\color{red}{\frac{2e}{3}}.$$

We have $$\sum_{m\geq1}\sum_{n\geq1}\frac{mn}{\left(n+m\right)!}=\sum_{m\geq1}\frac{1}{\left(m-1\right)!}\sum_{n\geq1}\frac{n}{\left(n-1\right)!}\int_{0}^{1}t^{m}\left(1-t\right)^{n-1}dt$$ $$=e\int_{0}^{1}t\left(2-t\right)dt=\color{red}{\frac{2}{3}e}$$ where the exchange of the series and the integral is justified by the dominated convergence theorem.

Using generating functions: let's start off with $$F(x, y) = \sum_{m=0}^\infty \sum_{n=0}^\infty \frac{x^m y^n} {(m+n)!}.$$ If we group the terms with $m+n = k$, we get $\frac{1}{k!} (x^k + x^{k-1} y + \cdots + y^k) = \frac{1}{k!} \cdot \frac{x^{k+1} - y^{k+1}}{x-y}$. Therefore, $$F(x,y) = \sum_{k=0}^\infty \frac{1}{k!} \cdot \frac{x^{k+1} - y^{k+1}}{x-y} = \frac{x e^x - y e^y}{x-y}.$$ Now, $$F_{xy}(x,y) = \sum_{m=1}^\infty \sum_{n=1}^\infty \frac{mn}{(m+n)!} x^{m-1} y^{n-1}.$$ Thus, the original sum is equal to $F_{xy}(1, 1)$. From here, it would be a straightforward but tedious calculation to find a closed-form formula for $F_{xy}$. That will give a fraction with $(x-y)^3$ in the denominator; however, it is not hard to see that $F_{xy}$ should be continuous at $(1,1)$, so you can calculate the value at $(1,1)$ as a limit of this quotient as $(x,y) \to (1,1)$. (For example, first substituting $x=1$ since $\{ (1,y) \mid y \ne 1 \}$ has cluster point at $(1,1)$, and then either using l'Hopital's rule three times or expanding a Taylor series about $y=1$ should work.) This will give the final answer $\frac{2}{3} e$.

For instance, in a Maxima session, I get (with a little hand editing of the transcript):

(%i1) diff(diff((x*exp(x)-y*exp(y))/(x-y),x),y);

(%i2) limit(limit(%o1,x,1),y,1);
          2 %e
(%o2)     ----
           3