I would use branching process (which is discrete-time random walk) to solve this problem.
According to your description, an individual "gives birth to" $2$ individuals with the probability $\frac{\lambda}{\lambda + \mu}$; an individual "gives birth to" $0$ individuals with the probability $\frac{\mu}{\lambda + \mu}$.
This gives the generating function $$G(s) = \frac{\mu}{\lambda + \mu} + \frac{\lambda}{\lambda + \mu}s^2$$.
Extinction probability is the smaller root of the function
$$s = \frac{\mu}{\lambda + \mu} + \frac{\lambda}{\lambda + \mu}s^2$$.
After calculation, the equation has two roots: $1$, $\frac{\mu}{\lambda} $.
Indeed, if $\frac{\mu}{\lambda} \geq 1$, i.e. $\mu \geq \lambda$, the population will go extinct within a finite time(𝑃 [ 𝑌𝑛=$0$ , 𝑛<∞].=$1$ )
A little bit more details. What is the probability that in a time interval $(0,t)$ exactly one person will have a success. Well first pick the person (say P1) that it will be call their child P2 and everyone else E, and condition on exactly when the success happens. Call that pdf of the first birth time of P1 $p_1(s)$.
$$P(\text{P1 has 1 birth & P2 and E have no births})$$
$$ = \int_0^t P(\text{P1 and P2 and E no birth after s}|\text{P1 has birth at time s})p_1(s)\mathrm{d}s$$
$$= e^{-\lambda t (n_0 - 1)}\int_{0}^t (e^{-\lambda (t-s)})^2 \lambda e^{-\lambda(s)} \mathrm{d}s$$
$$= e^{-\lambda t (n_0 - 1)}e^{-2 \lambda t} \left(e^{\lambda t}-1\right)$$
$$= e^{- \lambda n_0 t} \left(1-e^{-\lambda t}\right)$$
Then since we can make $n_0$ choices for the parent we get
$$P(\text{One birth in time interval (0,t)}) = n_0 e^{- \lambda t n_0} \left(1-e^{-\lambda t}\right)$$
Suppose that starting with $n_0$ the formula for $k$ births in the time interval is
$$P_t^{(n_0)}(k) = {n_0+k-1 \choose k} (e^{- \lambda t})^{n_0} (1 - e^{-\lambda t})^{k}$$
We wish to compute the formula for $P_t^{(n_0)}(k+1)$ by conditioning on when the first birth happens. It's a pretty elementary probability problem to figure out the PDF for the first of $n_0$ exponential RVs is
$$p(s) = \lambda n_0 e^{-\lambda n_0 s}$$
So we get
$$P_t^{(n_0)}(k+1) = \int_0^t\lambda n_0 \binom{k+n_0}{k} \left(1-e^{-\lambda (t-s)}\right)^k e^{-\lambda n_0 s-\lambda (n_0+1) (t-s)}\mathrm{d}s$$
$$= \frac{n_0 \binom{k+n_0}{k} \left(1-e^{-\lambda t}\right)^{k+1} e^{-\lambda n_0 t}}{k+1}$$
$$= \binom{(k+1)+n_0-1}{k+1} \left(1-e^{-\lambda t}\right)^{k+1} e^{-\lambda n_0 t}$$
Pretty easy to see that if we convert back to your notation using just the total number of people $n$ instead of the number gained $k$, that this gives the right thing.
Best Answer
As the comments mention, this is indeed a continuous time version of the Ehrenfest model and there are papers written considering this problem like
The model is a birth-death process, births when there are $i$ balls in box I are at rate $(N-i)\lambda$ and deaths at rate $i\lambda$. These formulas cover the special cases too. When box I is empty then there are no deaths, and births happen at rate $N\lambda$, while when box I is full there are no births and deaths occur at rate $N\lambda$.
The balance equations (writing $\pi_i$ for the probability of being in state $i$) are $$ \begin{align} (N\lambda)\pi_0 &= (N\lambda)\pi_1\\ (N\lambda)\pi_1 &= \lambda \pi_0 + (N-1)\lambda \pi_2 \\ (N\lambda)\pi_2 &= 2\lambda \pi_1 + (N-2)\lambda \pi_3 \\ & \vdots\\ (N\lambda)\pi_{N-1} &= (N-1)\lambda \pi_{N-2} + \lambda \pi_N \\ (N\lambda)\pi_N &= N\lambda \pi_{N-1} \\ \end{align}$$
which you then need to solve along with the constraint that $\sum_{i=0}^N \pi_i = 1$ to find the stationary probability distribution.