Solve a first order PDE derived from exponential generating function The PDE comes from a counting problem:
Suppose there is a bag of $m$ red, $n$ blue balls, and each time one randomly removes one ball until there are only balls of one color in the bag. What is the expected number of balls left?
Let $f(m,n)$ be the expectation. We have
\begin{equation*}
f(m,n)=\frac{m}{m+n}f(m-1,n)+\frac{n}{m+n}f(m,n-1),\quad f(m,0)=m,\; f(0,n)=n.
\end{equation*}
Multiply by $x^my^n/(m!n!)$ and sum over $m,n\geq 1$ to obtain the PDE:
\begin{equation*}
xF_x+yF_y-(x+y)F=xe^x+ye^y,\quad F(x,y)=\sum_{m,n\geq 0}f(m,n)\frac{x^my^n}{m!n!}.
\end{equation*}
My question is how to solve it?
Remark. Of course, if you know or can guess out the formula $f(m,n)=\frac{m}{n+1}+\frac{n}{m+1}$, the original problem can be solved easily by induction. Knowing the formula, you can also find
\begin{equation*}
F(x,y)=\left(\frac{x}{y}+\frac{y}{x}\right)e^{x+y}-\frac{x}{y}e^x-\frac{y}{x}e^y.
\end{equation*}
It does satisfy the PDE.

You have
$$x F_{x} + y F_{y} = (x + y) F + x e^{x} + y e^{y}$$
Using the method of characteristics in the parameterisation invariant form, we have
$$\frac{dx}{x} = \frac{dy}{y} = \frac{dF}{(x + y) F + x e^{x} + y e^{y}}$$
Solving across the first equality gives
$$\frac{x}{y} = C_{1}$$
Using this result, we can rewrite the last ratio as
\begin{align}
\frac{dF}{(x + y) F + x e^{x} + y e^{y}} &= \frac{dF}{y (x/y + 1) F + y (x/y) e^{y (x/y)} + y e^{y}} \\
&= \frac{dF}{y (C_{1} + 1) F + C_{1} y e^{C_{1} y} + y e^{y}}
\end{align}
and solving across the second equality gives
\begin{align}
\frac{dy}{y} &= \frac{dF}{y (C_{1} + 1) F + C_{1} y e^{C_{1} y} + y e^{y}} \\
\implies F' - (C_{1} + 1) F &= C_{1} e^{C_{1} y} + e^{y} \\
\implies (e^{-(C_{1} + 1) y} F)' &= C_{1} e^{- y} + e^{- C_{1} y} \\
\implies e^{-(C_{1} + 1) y} F &= - C_{1} e^{- y} - C_{1}^{-1} e^{- C_{1} y} + C_{2} \\
\implies F &= - C_{1} e^{C_{1} y} - C_{1}^{-1}  e^{y} + C_{2} e^{(C_{1} + 1) y}
\end{align}
Replacing the constant $C_{1} = x/y$ and noting $C_{2} = g(C_{1})$ for some arbitrary differentiable function $g$, we have
$$F = - \frac{x}{y} e^{x} - \frac{y}{x} e^{y} + g \left( \frac{x}{y} \right) e^{x + y}$$
which is the result you found, with $g(x/y) = x/y + y/x$. Note that $g$ is determined by the boundary conditions.
