Diffusion total flux with Finite Differences

Hi all,

I am trying to set a FD scheme to solve second Ficks law. I am particularly interested into getting the cumulative mass which has diffused out of my system (sphere). I am starting with the simplest situation (boundary concentration at r=R).

I have discretised the t and r domain and concentration profiles look right. I then calculate the total mass which has diffused out as 1-(C_average/C_initial).

When I compare my solution, it seems to be quite far from the analytical solution, given by Crank. I have increased the number of grid points but it does not make any difference, it stills under predict the total extracted mass with respect to the analytical solution.

I understand that a numerical method is never going to be as good as the analytical solution but still the error here is quite big..

Any help would be appreciated