% fd.m
% Solves a transient diffusion problem using the
% finite difference method.
global nt delx2

% set up the parameters
nt = 17
delx = 1/(nt-1)
delx2 = delx*delx
x(1) = 0.;
for i=2:nt-1
c0(i-1) = 0.;
x(i) = x(i-1)+delx;
end
x(nt) = x(nt-1)+delx

tbreak = [0 .001 .01 .1 1];
for kkk=1:4
tspan = [tbreak(kkk) tbreak(kkk+1)];
[t cc] = ode45('fd_rhs',tspan,c0);
nn=size(t);
c(kkk,:) = cc(nn(1),:);
c0(:) = c(kkk,:);
end

% add in the first and last point
for kkk=1:4
d(kkk,1) = 1;
for i=2:nt-1
d(kkk,i)=c(kkk,i-1);
end
d(kkk,nt) = 0;
end

plot(x,d(4,:),'g-*',x,d(3,:),'m-*',x,d(2,:),'b-*',x,d(1,:),'r-*')
xlabel('x')
ylabel('c')
legend('t = 1.0','t = 0.1','t = 0.01','t = 0.001')
title('Diffusion Problem Solved using Finite Difference Method')

******************************************************************

% fd_rhs.m
function ydot=fd_rhs(time,y)
global nt delx2

% add in the boundary conditions
c(1) = 1;
c(nt) = 0;

% transfer the y to the c
for i=1:nt-2
c(i+1) = y(i);
end

% calculate the right-hand side
for i=1:nt-2
dfd(i) = c(i) - 2*c(i+1) + c(i+2);
end
dfd = dfd/delx2;

% put result in the ydot
ydot = dfd';