# Solving diffusion equation using finite difference method

I am solving an 1-dimensional diffusion equation with Neumann boundary condition at outlet and constant concentration, C, at the inlet. In the end, I want to observe how the concentration diffuses over time along the x direction.

$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}$$

The following is the code written in MATLAB.

function main()
C0 = [500 480 460 440 420 400 380 360 340];
tspan = [0 10];
[t C]  = ode23s(@(t,s) fun(t,s), tspan , C0);
plot(t,C)

function dC = fun(t,C)
D  = 3000;
retain_rows = 2:8;
tri = -(D/(1/5^2))*full(gallery('tridiag',9,-1,2,-1)); % 5 is del x
tri = tri(retain_rows,:);
dC(1,1) = 0; % constant concentration at inlet
dC(retain_rows,1) = tri*C;
dC(9,1) = -D*((C(9)-C(8))/2); Neumann boundary condition at outlet
end
end


I have the solution obtained from ode solver. I'd like to know how to check for mole balance from the solution that is obtained. Also, could someone explain how the above code can be changed to apply constant flux at the inlet?

Any suggestions?