Can double generating functions be used to solve this? I was really interested in this question here
What is the expected number of red apples left when all the green apples are picked?

We have 4 green apples and 60 red apples. Each time we pick one out
  without replacement. Then what is the expected number of red apples
  left when all 4 green apples are picked?

I tried to use generating functions to get a closed form. I wrote a recurrence for the expected number of red apples like this (and I've verified that it is correct):
$$E(r,g) = \frac{r}{r+g}E(r-1,g) + \frac{g}{r+g}E(r,g-1)$$
and $E(0,g) = 0$ and $E(r,0) = r$.
So I did this
$$G(x) = \sum_{r=0}^{\infty} \sum_{g=0}^{\infty} E(r,g) x^r y^g$$
and tried expanding it out and shifting indices and taking derivatives but it was a huge mess and I had no confidence in what I was doing. Can generating functions be used here to get the closed form $E(r,g) = \frac{r}{g+1}$? I noticed this pattern by inspection and used induction to prove that the closed form was correct, but I wanted to know how to actually derive it using generating functions.
 A: [2016-08-13]:  A different way to proceed from the  PDE in (5), independent of other answers is stated  here. Coefficient extraction added at the end of this answer.

OPs recurrence relation is fine (if we additionally specify the range of validity) and we can derive a generating function from it.

We start by recalling the recurrence relation
  \begin{align*}
(r+g)E(r,g)&=rE(r-1,g)+gE(r,g-1)\qquad\qquad &r,g\geq 1\tag{1}\\
E(r,0)&=r &r\geq 0\\
E(0,g)&=0 &g\geq 0\\
\end{align*}

In the following we consider each term of the recurrence relation separately. 
It is also convenient to use the differential operator $D_x:= \frac{d}{dx}$ and $D_y:=\frac{d}{dy}$.

With the Ansatz
  \begin{align*}
G(x,y)=\sum_{r=0}^\infty\sum_{g=0}^\infty E(r,g)x^ry^g
\end{align*}
we obtain
  \begin{align*}
\sum_{r=1}^\infty&\sum_{g=1}^\infty(r+g)E(r,g)x^ry^g\tag{2}\\
&=(xD_x+yD_y)\sum_{r=1}^\infty\sum_{g=1}^\infty E(r,g)x^ry^g\tag{3}\\
&=(xD_x+yD_y)\left(G(x,y)-\sum_{g=0}^\infty E(0,g)y^g-\sum_{r=0}^\infty E(r,0)x^r+E(0,0)\right)\\
&=(xD_x+yD_y)\left(G(x,y)-\sum_{r=0}^\infty rx^r\right)\\
&=(xD_x+yD_y)\left(G(x,y)-(xD_x)\frac{1}{1-x}\right)\\
&=(xD_x+yD_y)G(x,y)-(xD_x)^2\frac{1}{1-x}\tag{4}\\
\end{align*}

Comment:


*

*In (2) we have to respect the validity of the recurrence relation and start with $r=g=1$ accordingly

*In (3) we use the operator $(xD_x)$ to transform a series $\sum_{r=1}^\infty a_r x^r$ to $\sum_{r=1}^\infty r a_r x^r$

*In (4) we use the linearity of the differential operator and since the geometric series is a function in $x$ the operator $yD_y$ transforms it to $0$.
Now we consider the the first summand of the RHS of (1).

\begin{align*}
\sum_{r=1}^\infty&\sum_{g=1}^\infty rE(r-1,g)x^ry^g\\
&=xD_x\sum_{r=1}^\infty\sum_{g=1}^\infty E(r-1,g)x^ry^g\\
&=xD_x\sum_{r=0}^\infty\sum_{g=1}^\infty E(r,g)x^{r+1}y^g\\
&=xD_xx\left(G(x,y)-\sum_{r=0}^\infty E(r,0)x^r\right)\\
&=xD_xx\left(G(x,y)-\sum_{r=0}^\infty r x^r\right)\\
&=xD_xxG(x,y)-xD_xx^2D_x\frac{1}{1-x}\\
\end{align*}

and  now the other summand of the RHS of (1).

\begin{align*}
\sum_{r=1}^\infty&\sum_{g=1}^\infty gE(r,g-1)x^ry^g\\
&=yD_y\sum_{r=1}^\infty\sum_{g=1}^\infty E(r,g-1)x^ry^g\\
&=yD_y\sum_{r=1}^\infty\sum_{g=0}^\infty E(r,g)x^ry^{g+1}\\
&=yD_yy\left(G(x,y)-\sum_{g=0}^\infty E(0,g)y^g\right)\\
&=yD_y yG(x,y)
\end{align*}

Putting all three expressions together we obtain

\begin{align*}
(xD_x+yD_y)&G(x,y)-(xD_x)^2\frac{1}{1-x}\\
&=xD_xxG(x,y)-xD_xx^2D_x\frac{1}{1-x}+yD_y yG(x,y)\\
(xD_x+yD_y&-xD_xx-yD_y y)G(x,y)\\
&=\left((xD_x)^2-xD_xx^2D_x\right)\frac{1}{1-x}\\
(xD_x(1-x)&+yD_y(1-y))G(x,y)=(xD_x)(1-x)(xD_x)\frac{1}{1-x}
\end{align*}
  and finally
  \begin{align*}
(xD_x(1-x)&+yD_y(1-y))G(x,y)=\frac{x}{(1-x)^2}\tag{5}
\end{align*}

We see $G(x,y)$ fulfills a linear differential equation in two variables. Please find   here a different way to proceed. But in fact we need not necessarily try to solve it, since we already know from  @MarkoRiedels answer:
\begin{align*}
E(r,g)=\frac{r}{g+1}
\end{align*}

So, we do the easy way and calculate $G(x,y)$ this way
\begin{align*}
G(x,y)&=\sum_{r=0}^\infty\sum_{g=0}^\infty E(r,g)x^ry^g\\
&=\sum_{r=0}^\infty\sum_{g=0}^\infty \frac{r}{g+1}x^ry^g\\
&=\left(\sum_{r=0}^\infty rx^r\right)\left(\sum_{g=0}^\infty \frac{1}{g+1}y^g\right)\\
&=xD_x\frac{1}{1-x}\cdot\frac{1}{y}\int\frac{1}{1-y}\,dy\\
&=-\frac{x}{(1-x)^2}\cdot\frac{1}{y}\left(\ln(1-y)+C\right)\tag{6}\\
\end{align*}

In order to determine the integration constant $C$ we evaluate $G(x,y)$ at $\left(\frac{1}{2},\frac{1}{2}\right)$. 

We obtain
  \begin{align*}
G\left(\frac{1}{2},\frac{1}{2}\right)
&=\sum_{r=0}^\infty r\left(\frac{1}{2}\right)^r\sum_{g=0}^\infty \frac{1}{g+1}\left(\frac{1}{2}\right)^g
=2\ln 4=4\ln 2
\end{align*}
  It follows together with (6)
  \begin{align*}
-4\left(\ln\frac{1}{2}+C\right)=4 \ln 2\qquad\longrightarrow\qquad C=0
\end{align*}
We finally conclude
  \begin{align*}
G(x,y)=-\frac{x}{(1-x)^2}\cdot\frac{\ln(1-y)}{y}
\end{align*}
  is a generating function for the sequence $\left(E(r,g)\right)_{r,g\geq 0}$ which fulfills the PDE (5).

[Update 2016-08-13]: According to OPs comments we extract the coefficient $E(r,g)$ from $G(x,y)$.
In order to do so, it's convenient to use the coefficient of operator $[x^r]$ to denote the coefficient of $x^r$ in a series $A(x)$.

We obtain for $r,g\geq 0$
  \begin{align*}
E(r,g)&=[x^ry^g]G(x,y)\\
&=[x^ry^g]\left(-\frac{x}{(1-x)^2}\cdot\frac{\ln(1-y)}{y}\right)\\
&=-[x^{r-1}y^{g+1}]\frac{\ln(1-y)}{(1-x)^2}\tag{1}\\
&=[x^{r-1}y^{g+1}]\sum_{k=0}^\infty\binom{-2}{k}(-x)^k\sum_{j=1}^\infty{\frac{1}{j}y^j}\tag{2}\\
&=\frac{1}{g+1}[x^{r-1}]\sum_{k=0}^\infty\binom{k+1}{1}x^k\tag{3}\\
&=\frac{r}{g+1}\tag{4}
\end{align*}

Comment:


*

*In (1) we use the rule $[x^p]x^qA(x)=[x^{p-q}]A(x)$.

*In (2) we apply the binomial series expansion and the use the series representation of $\ln(1-y)$.

*In (3) we select the coefficient of $y^{g+1}$ and use the binomial identity $\binom{-p}{q}=\binom{p+q-1}{p-1}(-1)^q$.

*In (4) we select the coefficient of $x^{r-1}$.
A: Suppose we have $G$ green apples  and $R$ red apples and ask about the
number of red apples left when all green apples have been picked.
In computing  a recurrence  for the expectation  $E(R, G)$  we note
that if the  sequence of samples ends in a $G$  the expected number of
red apples  left is  zero. If it  ends in  a $R$ it  is one  more than
$E(R-1, G).$ We obtain the recurrence
$$E(R, G) = {R+G\choose R}^{-1} {R-1+G\choose R-1}
(1 + E(R-1, G))$$
where $E(0, G) = 0.$ Now
$${R+G\choose R}^{-1} {R-1+G\choose R-1} = 
\frac{(R-1+G)! R! G!}{(R+G)! (R-1)! G!}
= \frac{R}{R+G}$$
so this is
$$E(R, G) = \frac{R}{R+G}
(1 + E(R-1, G)).$$
We may drop the $G$ as an index variable as it never changes. We get
$$R E(R) + G E(R) = R + R E(R-1).$$
Introducing the generating function
$$F(z) = \sum_{R\ge 0} E(R) z^R$$
and multiplying by $z^R$ and summing over $R\ge 1$ we get
$$\sum_{R\ge 1} R z^R [z^R] F(z)
+ \sum_{R\ge 1} G z^R [z^R] F(z)
= \sum_{R\ge 1} R z^R
+ \sum_{R\ge 1} R z^R [z^{R-1}] F(z).$$
This is
$$z F'(z) + G F(z) = \frac{z}{(1-z)^2} +
\sum_{R\ge 1} R z^R [z^{R-1}] F(z).$$
Solving the remaining term we get
$$\sum_{R\ge 2} R z^R [z^{R-1}] F(z)
= z^2 \sum_{R\ge 2} R z^{R-2} [z^{R-1}] F(z)
\\ = z^2 \sum_{R\ge 2} (R-1) z^{R-2} [z^{R-1}] F(z)
+ z^2 \sum_{R\ge 2} z^{R-2} [z^{R-1}] F(z)
\\ = z^2 F'(z) + z F(z).$$
This ODE is linear and we get
$$F'(z) + \frac{G-z}{z(1-z)} F(z) = \frac{1}{(1-z)^3}.$$
The integrating factor is
$$\exp((1-G)\log(z-1) + G\log(z))
= z^G (z-1)^{1-G}$$
and we obtain
$$F(z) = z^{-G} (z-1)^{G-1}
\left(\kappa - \int \frac{1}{(z-1)^3} z^G (z-1)^{1-G} dz\right)
\\ = z^{-G} (z-1)^{G-1}
\left(\kappa + \frac{1}{G+1} \frac{z^{G+1}}{(z-1)^{G+1}}\right)
\\ = \kappa \frac{1}{z} \left(1-\frac{1}{z}\right)^{G-1}
+ \frac{1}{G+1} \frac{z}{(z-1)^2}.$$
It follows even without  determining $\kappa$ (the polynomial in front
contains only negative powers of $z$) that for $R\ge 0$
$$\bbox[5px,border:2px solid #00A000]{
[z^R] F(z) = \frac{1}{G+1} [z^R] \frac{z}{(1-z)^2}
= \frac{R}{G+1}.}$$
