Partial Differential Equation to evaluate double series I was trying to compute a closed form for the double series
$$f(x,y,z)=\sum_{m,n=1}^\infty \frac{x^m y^n z^{m+n}}{m!(m+n)!}$$
and I noticed the fact that
$$\frac{\partial^2f}{\partial x\partial z}-f=e^y-1$$
I know how to solve some somewhat advanced differential equations for functions of one variable, but I know nothing about how to solve partial differential equations. Can someone show me how to solve this one?
I did have one little insight about this equation. Since no derivatives are taken with respect to $y$ on the LHS of the differential equation, I can just treat $e^y-1$ as if it was a constant and instead solve the differential equation as if $f$ were a function $f(x,z)$ of two variables only. But since I don't know how to get started, this doesn't help me.
 A: I do not know if a closed form could be found for 
$$S_{m,n}=\sum_{m=1}^\infty \sum_{n=1}^\infty \frac{x^m\, y^n\, z^{m+n}}{m!\,(m+n)!}=\sum_{m=1}^\infty \sum_{n=1}^\infty \frac{(xz)^m \,(yz)^n }{m!\,(m+n)!}=\sum_{m=1}^\infty \sum_{n=1}^\infty \frac{a^m\, b^n }{m!\,(m+n)!}$$ but, at least, the inner sum can be computed since 
$$\sum_{n=1}^\infty \frac{a^m\, b^n }{m!\,(m+n)!}=\frac{a^m}{m!}\sum_{n=1}^\infty \frac{ b^n }{(m+n)!}=\frac{1}{m!}\left(\frac ab\right)^m\, e^b\left(1-\frac{\Gamma (m+1,b)}{\Gamma (m+1)} \right)$$ 
A: $$f(x,y,z)=\sum_{m=1}^\infty\sum_{n=1}^\infty  \frac{x^m y^n z^{m+n}}{m!(m+n)!}$$
$$\frac{\partial f}{\partial z}=\sum_{m=1}^\infty\sum_{n=1}^\infty  \frac{x^m y^n z^{m+n-1}}{m!(m+n-1)!}$$
Let $n=\nu+1\quad;\quad \frac{\partial f}{\partial z}=\sum_{m=1}^\infty\sum_{\nu=0}^\infty  \frac{x^m y^{\nu+1} z^{m+\nu}}{m!(m+\nu)!}$
$$\frac{\partial f}{\partial z}=y\sum_{m=1}^\infty\sum_{\nu=1}^\infty  \frac{x^m y^{\nu} z^{m+\nu}}{m!(m+\nu)!} +y\sum_{m=1}^\infty  \frac{x^m  z^{m}}{m!m!}$$
$\sum_{m=0}^\infty  \frac{(xz)^m}{(m!)^2}=I_0(2\sqrt{xz}) \quad\to\quad \sum_{m=1}^\infty  \frac{(xz)^m}{(m!)^2}=I_0(2\sqrt{xz})-1 $
$I_0$ is the modified Bessel function of first kind and order $0$.
$$\frac{\partial f}{\partial z}=yf(x,y,z)+y\left(I_0(2\sqrt{xz})-1\right)$$
This is a first order linear differential equation with respect to the variable $z$ and in which $x,y$ are parameters. Solving it is straightforward with usual method. The result is :
$$f(x,y,z)=1+ye^{yz}\left(\int_0^z e^{-y\xi}I_0(2\sqrt{x\xi})d\xi+F(x,y)\right)$$
The arbitrary $F$ is determined according to the particular value :
$\quad f(x,y,0)=0=1+yF(x,y)\quad$ which leads to :
$$f(x,y,z)=1-e^{yz}+ye^{yz}\int_0^z e^{-y\xi}I_0(2\sqrt{x\xi})d\xi$$
As far as I know, there is no closed form for the integral. As a consequence, the final result involves a function defined by an integral.
$$f(x,y,z)=\sum_{m=1}^\infty\sum_{n=1}^\infty  \frac{x^m y^n z^{m+n}}{m!(m+n)!}=1-e^{yz}+ye^{yz}\int_0^z e^{-y\xi}I_0(2\sqrt{x\xi})d\xi$$
NOTE :
Obviously the result is valid for $y=0$ or $z=0$. It is also valid for $x=0$ since $I_0(0)=1$ then $f(0,y,z)=1-e^{yz}+ye^{yz}\int_0^z e^{-y\xi}d\xi=1-e^{yz}+ye^{yz}\left(-\frac{e^{-yz}-1}{y}\right)=0$
