When one considers the 1-D diffusion equation in cartesian coordinates
$$\frac{\partial u}{\partial t}=\chi\frac{\partial^2u}{\partial x^2},$$
one finds that the amplification factor for the Crank Nicolson scheme (with central differences in the spatial derivatives) is
$$A=\frac{1-2F\sin^2(k\Delta x/2)}{1+2F\sin^2(k\Delta x/2)},$$
where $F=\chi\Delta t/(\Delta x)^2$. Meaning that this method is unconditionally stable, i.e., the method is stable whatever $F$ is.
Now, let us consider the equation in cylindrical coordinates with a source term:
$$\frac{\partial u}{\partial t}=\chi\frac{1}{r}\frac{\partial}{\partial r}\left[r\frac{\partial u}{\partial r}\right]+S_{ext}$$
with boundary conditions $\partial_rT(r=0)=0$ and $T(r=1)=T_a$.
When $S_{ext}=0$, I get that the amplification factor is given by
$$A=\frac{1-\frac{F}{2}\left[4\sin^2\left(\frac{k\Delta r}{2}\right)-\frac{\sin(k\Delta r)}{q}\right]}{1+\frac{F}{2}\left[4\sin^2\left(\frac{k\Delta r}{2}\right)-\frac{\sin(k\Delta r)}{q}\right]}$$
where $q$ is the cell number. I am not sure if this is the correct result, namely it is strange to me the fact that it depends on $q$. Indeed, for small $q$ it is easy to see that $A$ may become larger than 1 for some of the wavenumbers $k$ and independently of $F$ .
Questions:
1. Does this mean that the method may become unstable?
2. When onde considers $S_{ext}\neq0$ the same result for $A$ applies?
3. And what about the accuracy in the $S_{ext}\neq0$ case? Because I have found that, with increasing $S_{ext}$, the deviation from the theoretical result increases a lot before the steady-state is reached.
The Von Neumann stability analysis only works for linear PDEs with constant coefficients. It is not the case for your problem because of the $r$ and $1/r$ terms.