$$ \sum_{n=0}^{\infty}\left(\frac{x^n}{n!}\sum_{r=0}^{n-1}\left(\frac{y^r}{r!}\right)\right) $$
This is a double summation I need to evaluate(not for a homework problem, but for a probability question involving gaming I found). I can't find any idea how to do this and wolfram alpha can't calculate this with extended computer time either, though this series is obviously less than e^(x+y) and monotonous, so it converges.
I will show you how to evaluate the special case for $x=y$ in terms of a modified Bessel function of the first kind. I will also give you a functional equation for your series in terms of the exponential function and a modified Bessel function of the first kind.
If we denote your function as:
$$f(x,y)=\sum_{n=0}^{\infty}\left(\frac{x^n}{n!}\sum_{r=0}^{n-1}\left(\frac{y^r}{r!}\right)\right)$$
And the modified Bessel function of the first kind with index zero as:
$$I_0(x)=\sum_{n=0}^\infty\frac{x^{2n}}{4^nn!^2}$$
Then we have that,
$$f(x,x)=\frac{e^{2x}-I_0(2x)}{2}$$
And in addition,
$$f(x,y)+f(y,x)=e^{x+y}-I_0(2\sqrt{xy})$$
To prove this we will define a slightly simpler function $h$ for convenience: $$h(x,y)=\sum_{n=0}^\infty\frac{x^n}{n!}\sum_{k=0}^n\frac{y^k}{k!}$$
Then we have that,
$$h(x,y)=\sum_{n=0}^\infty\frac{x^n}{n!}\sum_{k=0}^n\frac{y^k}{k!}=\sum_{n=0}^\infty\frac{x^n}{n!}\sum_{a+b=n}\frac{y^b}{b!}=\sum_{n=0}^\infty\sum_{a+b=n}\frac{x^{a+b}y^b}{(a+b)!b!}=\sum_{a=0}^\infty\sum_{b=0}^\infty\frac{x^{a+b}y^b}{(a+b)!b!}$$
$$=\sum_{a=0}^\infty\sum_{b=0}^\infty\frac{y^bx^{a+b}}{b!(a+b)!}=\sum_{b=0}^\infty \frac{y^b}{b!}(\frac{x^{b}}{b!}+\frac{x^{1+b}}{(1+b)!}+\frac{x^{2+b}}{(2+b)!}+\frac{x^{3+b}}{(3+b)!}+\dots)$$
$$h(x,y)=\sum_{b=0}^\infty \frac{y^b}{b!}(\frac{x^{b}}{(b)!}+\frac{x^{1+b}}{(1+b)!}+\frac{x^{2+b}}{(2+b)!}+\frac{x^{3+b}}{(3+b)!}+\dots)$$ $$h(y,x)=\sum_{b=0}^\infty\frac{y^b}{b!}(\frac{x^0}{0!}+\frac{x^1}{1!}+\frac{x^2}{2!}+...\frac{x^{b-1}}{(b-1)!}+\frac{x^b}{b!}+\dots)$$
$$h(x,y)+h(y,x)=\sum_{b=0}^\infty \frac{y^b}{b!}(\frac{x^b}{b!}+e^x)$$ $$h(x,y)+h(y,x)=e^{x+y}+I_0(2\sqrt{xy})$$
$$h(x,y)=\sum_{n=0}^{\infty}\left(\frac{x^n}{n!}\sum_{r=0}^{n-1}\left(\frac{y^r}{r!}\right)\right)+\sum_{n=0}^\infty\frac{(xy)^n}{n!^2}=f(x,y)+I_(2\sqrt{xy})$$
Thus by simplification,
$$f(x,y)+f(y,x)=e^{x+y}-I_0(2\sqrt{xy})$$
And the case for $x=y$ follows by simple algebraic manipulation.
The only thing I was able to find is the closed form expression of the inner summation. The result is given by $e^y \Gamma(n, y) / \Gamma(n)$ where $\Gamma(n, y)$ is the incomplete Gamma function. As you noticed, this is smaller than $e^y$ and the complete series is less than $e^{x+y}$.
I am not sure that a closed form could be found.
Happy New Year !
