$$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)!}$$ Dejemos que $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$ es la función de Bessel modificada de primer tipo y orden $0$ . $$\frac{\partial f}{\partial z}=yf(x,y,z)+y\left(I_0(2\sqrt{xz})-1\right)$$ Se trata de una ecuación diferencial lineal de primer orden con respecto a la variable $z$ y en el que $x,y$ son parámetros. Resolverlo es sencillo con el método habitual. El resultado es : $$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)$$ La arbitrariedad $F$ se determina en función del valor particular : $\quad f(x,y,0)=0=1+yF(x,y)\quad$ que lleva a : $$f(x,y,z)=1-e^{yz}+ye^{yz}\int_0^z e^{-y\xi}I_0(2\sqrt{x\xi})d\xi$$ Que yo sepa, no hay una forma cerrada para la integral. En consecuencia, el resultado final implica una función definida por una 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$$ NOTA :
Obviamente, el resultado es válido para $y=0$ o $z=0$ . También es válido para $x=0$ desde $I_0(0)=1$ entonces $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$