Numerical solution of non-constant coefficient diffusion equation via finite-difference method I would like to solve numerically a diffusion equation:
$$\frac{\partial c}{\partial t} = \frac{\partial}{\partial x}\left(k\frac{\partial c}{\partial x}\right),$$
given some boundary conditions and initial value.
My idea is to solve it via discretization of the spatial variable using finite-difference method and then solve the resulting system of ODEs with Matlab, that is, use the Method of Lines approach. However, as the term $k=k(x)$ then I do not know how to discretize the right hand side of the equation. Could you suggest any idea?
Thanks!
 A: The straight forward way would be to approximate the right hand side by
$$
\frac1h\Bigl(k(x+h)\frac{c(x+2h)-c(x+h)}{h}-k(x)\frac{c(x+h)-c(x)}{h}\Bigr).
$$
A: If $k$ is continuously differentiable, then there is a simple approach:
$$\frac{\partial}{\partial x} \left ( k \frac{\partial c}{\partial x} \right ) = k'(x) \frac{\partial c}{\partial x} + k(x) \frac{\partial^2 c}{\partial x^2}.$$
Now you can discretize these two spatial derivatives as you see fit.
More generally, it is better to think about fluxes. (Indeed it is essential to think about fluxes in hyperbolic problems like the transport and wave equations, even though it is not so crucial in parabolic problems like the diffusion equation.) Consider a point $x$ on the line with neighbors $x-\Delta x$ and $x+\Delta x$. The flux balance at this point over a time step of length $\Delta t$ should look like this:
\begin{align}
x-\Delta x \to x & : k(x-\Delta x/2) \Delta t c(x-\Delta x) \\
x \to x-\Delta x & : k(x-\Delta x/2) \Delta t c(x) \\
x \to x+\Delta x & : k(x+\Delta x/2) \Delta t c(x) \\
x+\Delta x \to x & : k(x+\Delta x/2) \Delta t c(x+\Delta x)\end{align}
(plus corrections of higher order, of course). 
This would give the ODE form
$$\frac{\partial c}{\partial t}(t,x)=k(x-\Delta x/2)c(t,x-\Delta x)+k(x+\Delta x/2)c(t,x+\Delta x) \\
-c(t,x)(k(x-\Delta x/2)+k(x+\Delta x/2)).$$
If $k$ is continuous then this will work, though it could be slow if $k$ is especially irregular. If $k$ is not even continuous, then you must be even more careful than this, and will need to write down a weak formulation of the problem in order to write down a sensible method.
