% solvediff.m
% Numerically solve the one dimensional diffusion equation
% Numerical method: central differential method
% 06/26/2009 created by Xiaozhou Zhou
X=linspace(0,1,101); % meshgrid
imax=length(X); % number of mesh elements
dx=X(2:end)-X(1:end-1); % vector of element width
dxmin=min(dx); % minimum width
tttl=3;
dt=dxmin^2/4; % time step
nt=floor(tttl/dt); % number of time steps
T=0:dt:nt*dt;
C0 = zeros(1,imax); % initial condition
C=C0;
iout=1;
nout=10000;
Cs(iout,:)=C;
for it=1:nt
C1=cds(C,imax-1,dx,dt);
C=C1;
% boundary conditions
C(1)=1;
C(imax)=0;
% output
if mod(it,nout)==0
iout=iout+1;
Cs(iout,:)=C1;
end
end
% cds.m
% central differential method
function fun=cds(C,m,dx,dt)
fun=C;
for j=2:m
fun(j)=C(j)+2*dt/(dx(j-1)+dx(j))*((C(j+1)-C(j))/dx(j)-(C(j)-C(j-1))/dx(j-1));
end