I'm trying to approximate an axial-symmetric diffusion equation in cylindric coordinates in the domain $0 \leq r \leq R, 0 \leq z \leq l$ using Crank-Nicolson scheme and facing a sort of wierd and wild behaviour at r=0.
The scheme is as follows:
$$ \frac{u^{k+1}-u^k}{\tau} = \frac{D}{2}((\frac{\frac{i+\frac{1}{2}}{i}u_{i+1}+\frac{i-\frac{1}{2}}{i}u_{i-1}-2iu_{ij}}{ih^2})^k + (\frac{\frac{i+\frac{1}{2}}{i}u_{i+1}+\frac{i-\frac{1}{2}}{i}u_{i-1}-2iu_{ij}}{ih^2})^{k+1}) + \frac{D}{2}((\frac{u_{j+1}+u_{j-1}-2u_{ij}}{h^2})^k + (\frac{u_{j+1}+u_{j-1}-2u_{ij}}{h^2})^{k+1}) $$
At r=0:
$$ \Delta = \Delta_r + \Delta_z $$
$$ \Delta_r = \frac{4}{h^2}(u_1 - u_0) $$
For the uniform and exponential distributions at t=0
i am getting negative values for $u_{0j}^k$ at timesteps k=1,2,3 ....
Is it anyhow avoidable?