I am trying to use a explicit FDM for transient 1d conditions with linear elements for specific diffusion equation: ds/dt=D(x)*d2s/dx2+R(s) the problem is that I am using different, nonlinear functions describing diffusion constant D(x) and reaction rate R(s). The diffusion parameter is dependent on position along x and in general case is a function with non-continuous dD(x)/dx at some points. I am using average value of diffusion parameter for each linear element but the solution seems to be not right. I think that is the problem, because when I use constant D everything looks all right. Especially plot of d2s/dx2 is very wired. What I am doing wrong? Thanks.