Solution (source code)

= Solution

For $f_n=n^2$, the birth increment is $2n+1$, while
$$
f_n-f_{n-2}=n^2-(n-2)^2=4n-4.
$$
With $D=2d$, the generator identity gives
$$
\boxed{
\frac d{dt}\langle n^2\rangle
=b(2\mu+1)-4d(\langle n^2\rangle-\mu)-2dp_1}.
$$
At stationarity $\mu=b/(2d)$, and solving this equation for the <variance> yields
$$
\sigma^2=\langle n^2\rangle-\mu^2
=\frac{3b}{4d}-\frac{p_1}{2}.
$$
For $b,d>0$ the chain is irreducible on the nonnegative integers, so its stationary distribution has $0<p_1<1$. Therefore
$$
\boxed{\frac{3b}{4d}-\frac12<\sigma^2<\frac{3b}{4d}}.
$$