Let $X_t$ be a simple birth process (with birth rate $\lambda_n = n\lambda$) such that $X_0 = 1$. Set $\mu_t = E[X_t]$ and $\nu_t = E[X_t^2]$. Find differential equations for $\mu_t$ and $\nu_t$ and hence calculate the variance of $X_t$.
I've tried using the Kolmogorov forward differential equation to get $\mu_t$ which gave me $e^{-\lambda t}$ but I'm not sure how to calculate $\nu_t$.
The differential equations (with initial conditions) for this Yule process with linear birth rates are \begin{align} P_1'(t) &= -\lambda P_1(t),\quad P_1(0)=1\tag 1\\ P_n'(t) &= -n\lambda P_n(t)+(n-1)\lambda P_{n-1}(t),\quad P_n(0)=0.\\ \end{align} Solving $(1)$ yields $P_1(t) = e^{-\lambda t}$. Assume now that $P_n(t) = e^{-\lambda t}(1-e^{-\lambda t})^{n-1}$. Then \begin{align} P_{n+1}'(t) &= -(n+1)\lambda P_{n+1}(t) + n\lambda P_n(t)\\ &= -(n+1)\lambda P_{n+1}(t) + n\lambda e^{-\lambda t}(1-e^{-\lambda t})^{n-1}, \end{align} so $$ P_{n+1}'(t) +(n+1)\lambda P_{n+1}(t) = n\lambda e^{-\lambda t}(1-e^{-\lambda t})^{n-1}.\tag2 $$ This is a first-order linear differential equation, with integrating factor $$ \mu(t) = e^{\int (n+1)\lambda \ \mathsf dt} = e^{(n+1)\lambda t}. $$ Multiplying $(2)$ by $\mu(t)$, we have $$ \frac{\mathsf d}{\mathsf dt}[P_{n+1}(t)e^{(n+1)\lambda t}] = n\lambda e^{n\lambda t}(1-e^{-\lambda t})^{n-1}. $$ Integrating and multiplying by $1/\mu(t)$ yields $$ P_{n+1}(t) = e^{-\lambda t}(1-e^{-\lambda t})^n,\tag3 $$ so by induction $(3)$ holds for all $n\geqslant 1$. The first moment of $X(t)$ is given by \begin{align} \mathbb E[X(t)] &= \sum_{n=1}^\infty n P_n(t)\\ &= \sum_{n=1}^\infty ne^{-\lambda t}(1-e^{-\lambda t})^{n-1}\\ &= e^{\lambda t}, \end{align} and the second moment by \begin{align} \mathbb E[X(t)^2] &= \sum_{n=1}^\infty n^2 P_n(t)\\ &= \sum_{n=1}^\infty n^2e^{-\lambda t}(1-e^{-\lambda t})^{n-1}\\ &= e^{\lambda t} \left(2 e^{\lambda t}-1\right), \end{align} so the variance is \begin{align} \mathsf{Var}(X(t) &= E[X(t)^2] - E[X(t)]^2\\ &= e^{\lambda t} \left(2 e^{\lambda t}-1\right) - (e^{\lambda t})^2\\ &= e^{\lambda t} \left(e^{\lambda t}-1\right). \end{align}