Crank-Nicolson solution of parabolic PDE with Newumann boundary conditions

I am trying to solve the non-linear parabolic PDE in $c(t,r)$
$$c_t=frac{1}{r}(rDc_r-alpha r^2 c)_r$$
with initial condition $c(0,r)=f(r)$
and boundary conditions $c_r(t,r_1)=alpha r_1c_1/D$ and $c_r(t,r_n)=alpha r_nc_n/D$; in other words, there is no flux across the boundaries.

This the Lamm equation, a 1-dimensional problem in cylindrical coordinates. In the simplest case, $D$ and $alpha$ are constants.

I have set up a uniformly-spaced $r$ grid with $n$ points, so $r_j=r_1+(j-1)Delta r$, with $1<=j<=J$. I’m using the Crank-Nicolson method with forward difference for $c_t$, and central differences for $c_{rr}$ and $c_r$. At the boundaries, I have added "ghost points" $c_0$ and $c_{J+1}$.

The difference equation resulting is, for $j=1,ldots, J$
$$frac{c_{n+1,j}-c_{n,j}}{Delta t}=
frac{D}{2}left[frac{c_{n,j+1}-2c_{n,j}+c_{n,j-1}}{(Delta r)^2}+frac{c_{n+1,j+1}-2c_{n+1,j}+c_{n+1,j-1}} {(Delta r)^2}right]
+left(frac{D}{ 2r_j}- frac{alpha r_j}{ 2}right) left[frac{c_{n,j+1}-c_{n,j-1}}{ 2Delta r}
+frac{c_{n+1,j+1}-c_{n+1,j-1}}{ 2Delta r}right]
-alphaleft[c_{n,j}+c_{n+1,j}right]$$
where now $c_{n,j}$ is the value of $c$ at the radial point $j$ and the $n$-th time step. Values for the "ghost points" are obtained from the above equation applied at $j=1$ and $j=J$ and the boundary conditions. They are then substituted into the difference equation.

This formulation leads to the matrix equation ${bf Fc}^{n+1}={bf Gc}^n$, where ${bf F}$ and ${bf G}$ are tridiagonal, and ${bf c}^k$ is a column vector of $c_{k,j}$ at time step $k$.

Solving for ${bf c}^{n+1}$ gives a result that does not conserve mass (i e, $c$). Since the boundaries have 0 flux, $int_{r_1}^{r_J} r dr$ should remain constant, but it does not.

Can anyone advise me if there is something I have done incorrectly, or where to look for an error in the computations? I did not include the elements of ${bf F}$ and ${bf G}$, but I could make those available.