I want to implement a finite difference code of this simple diffusion equation with space-dependent diffusivity:
$$\partial_{t}u =D\partial_{x}^{2}u+\partial_{x}D\cdot\partial_{x}u$$
I go for a fully implicit method for the diffusion part, and explicit for the first-order-derivatives term:
$$\frac{(\Delta x)^{2}}{\Delta t}(u_{n+1,j}-u_{n,j})=D(x_{j})[u_{n+1,j+1}-2u_{n+1,j}+u_{n+1,j-1}]\\+[D_{11}(x_{j+1})-D(x_{j})](u_{n,j+1}-u_{n,j})$$
where $n$ is the time index and $j$ runs over the grid points, $\Delta x$ is the grid spacing and $\Delta t$ is the time step.
Using a time-step halving technique to control the error (that is, essentially, I calculate $u_{n+1}$ with $\Delta t$ and $\Delta t/2$ and try to minimize the difference) I am sure that this scheme will be 2nd-order accurate in time and stable [see Richardson extrapolation].
I use as initial condition a bell-function with $u=0$ at both boundaries ($x=0,1$).
The code which results from this, does not seem to work, in fact it gives a strange shift toward the right boundary, as if something was wrong in the spacing or the velocity at which diffusion is happening? Can this be a problem of the scheme I am using?