Let's consider a process in which individuals die and give birth with rate one. This process is known as a continuous time critical branching process. The probability to observe \(n\) individual at time \(t\) obeys the following difference equation: $$ P(n, t+dt) = (1-2 n dt )P(n, t) + dt\left[ (n-1) P(n-1,t) + (n+1)P(n+1,t)\right] $$ where the first term corresponds to the case when nothing is happening in the time interval \(dt\), the second corresponds to division/birth and the third to the death of one individual with rate. Taking the limit \(dt\to 0\), we immediately get $$ \partial_t P(n, t) = -2nP(n, t) + (n-1)P(n-1,t) + (n+1)P(n+1,t) $$ To solve this equation, we define the generating function \(\hat{p}(z,t) = \sum_n z^n P(n,t)\), multiply the equation by \(z^n\) and sum over \(n\): $$ \partial_t \hat{p}(z, t) = -2z \partial_z \sum_{n=0} z^{n} P(n, t) + z^2 \partial_z \sum_{n=0} z^{n-1}P(n-1,t) + \partial_z \sum_{n=1} n z^{n+1}P(n+1,t) = \left[-2z + z^2 + 1\right]\partial_z \hat{p}(z,t) $$ Instead of the difference equation, we are now looking at partial differential equation for the generating function.
We will consider here the initial condition \(P(n,0)=\delta_{n0}\), i.e., there is exactly one individual at time \(t=0\), and \(\hat{p}(z,t) = z\) and solve this partial differential equation using the methods of characteristics. The method of characteristics determines families of curves in the space of \((z,t)\) along which the function is constant. Parameterizing the characteristic with \(t\), we have \(\hat{p}(z(t),t)\) and consequently. $$ 0 = \frac{d\hat{p}(z(t), t)}{dt} = \partial_z \hat{p}(z,t) \frac{dz}{dt} + \partial_t \hat{p}(z,t) = \partial_z \hat{p}(z,t) \left[\frac{dz}{dt} -(b+d)z + b z^2 + d \right] $$ Substituting \(z(t) = 1 - \chi(t)\), this equation simplifies to $$ \frac{d \chi(t)}{dt} = \chi(t)^2 $$ which has a finite time singularity. $$ \chi(t) = \frac{1}{C - t} $$ with initial condition \(\chi(0) = C^{-1}\). By construction, the generating function \(\hat{p}(z,t)\) is constant along the characteristic curve such that \(\hat{p}(z(t),t) = \hat{p}(z(0),0)) = z(0) = 1-C^{-1}\). The problem is hence reduced to finding the integration constant \(C\) such that the characteristic passes through the point \((z,t)\). Solving $$ z = 1 - \frac{1}{C(z,t)-t} $$ for \(C(z,t)\) yields $$ % (1-z)^{-1} = C(z,t)-t C(z,t) = (1-z)^{-1} + t $$ This in turn results in the generating function $$ \hat{p}(z,t) = 1 - \frac{1}{(1-z)^{-1} + t} $$ As a sanity check, we can evaluate the first few moments $$ 1 = \sum_n P(n,t) = \hat{p}(1,t) = 1 $$ $$ \langle n\rangle = \sum_n n P(n,t) = \partial_z \hat{p}(z,t)|_{z=1} = 1 $$ and the survival probability $$ 1-P(0,t) = 1-\hat{p}(0,t) = \frac{1}{1+t} $$
Solution via the first step equation
While the way we derived the generating function above was straightforward, it did require solving a (simple) PDE. A more direct way to the solution is via a first step equation. Splitting the time interval \([0,t]\) in \([0,\Delta t]\) and \([\Delta t, t]\), and assuming we start with a single individual at time \(t=0\), we have $$ P(n,t) = (1-2\Delta t)P(n,t-\Delta t) + \Delta t \left[\delta(n=0) + \sum_{n'=0}^n P(n',t)P(n-n',t)\right] $$ The first term corresponds to the case where nothing happens in the interval \([0,\Delta t]\), the third term corresponds to the death of the one and only individual (and pins the distribution at \(n=0\)), while the last term corresponds to a birth event and the sum runs over the different ways in which \(n\) descendants can result from two seeds.
This equation is readily turned into an ODE $$ \frac{dP(n,t)}{dt} = -2 P(n,t) + \delta(n=0) + \sum_{n'=0}^n P(n',t)P(n-n',t) $$ and the convolution can be disentangled by Laplace-transformation, that is considering the generating function \(\hat{p}(\lambda,t)\ = \sum_n \lambda^n P(n,t)\). The Laplace transform turns the convolution into a product: $$ \frac{d\hat{p}(\lambda,t)}{dt} = -2 \hat{p}(\lambda,t) + 1 + \sum_{n}\sum_{n'=0}^n \lambda^{n'}P(n',t)\lambda^{n-n'}P(n-n',t) = 2\hat{p}(\lambda,t) + 1 + \hat{p}(\lambda,t)^2 $$ (Note that \(\sum_{n=0}^\infty)\sum_{n'=0}^{n} = \sum_{n'=0}^\infty)\sum_{n=n'}^{\infty} \)). This can be further simplified by substituting \(\hat{p}(\lambda,t) = 1-\phi(\lambda,t)\): $$ \frac{d\phi(\lambda,t)}{dt} = 2 (1-\phi) - 1 - (1-\phi))^2 = -\phi^2 $$ with solution \(\phi(\lambda,t) = \frac{1}{(1-\lambda)^{-1}+t}) \), where the integration constant was used to fulfill the initial condition \(\phi(\lambda,t=0)=1-\lambda\).