 % This file solves the diffusion equation with Dirichlet BC's
 %using Crank Nicholson method for an exercise
clear

L = 1; %length of the domain
D = 1; % the diffusion coefficient
u0 = 0; % boundary value at x=0
uL = 0 ; % boundary value at x=L
 N = 20; %N+1 = number of grid points, N-1 = number of interior grid points
 % there are N-1 unknowns
 dx = L/(N); % grid size
 
 bc = [u0;zeros(N-3,1);uL]; % boundary contribution to difference operator
 sc =  2*ones(N-1,1) ; % diagonal of difference operator
 scu = [ones(N-2,1)]; %upper diagonal for difference operator 
 scl = [ones(N-2,1)]; %lower diagonal for difference operator
 
 X = [1:N-1]'*dx;  % interior grid points

 %
 
%pick a time step:
  
M = 40;  % number of time steps
 
dt = dx/5; % time step  %this is much bigger than allowed for forward euler
 
scal = D*dt/dx^2/2;  %The factor of 1/2 is because of the CN step split

  %V = exp(-(X-dx*N/2).^2./(2*dx));  % initial data is gaussian
  V = 2*sin(pi*X);  % initial data for exercise
kstep = 1;
t=0;

% uses Crank Nicolson to solve the diffusion equation

% set up the implicit matrix
 
Am=diag(1+ scal*sc)-scal*diag(scu,1) -scal*diag(scl,-1);
% The algorithm

for n=2:M
for k = 1:kstep
  
V_t = scal * ([scu;0].*[V(2:N-1 );0] - sc.*V + [0;scl].*[0;V(1:N-2 )]+2*bc);
V = Am\(V+V_t);
 
t=t + dt;
end

% exact solution
Vexact = 2*exp(-pi^2*t)*sin(pi*X);
figure(1)
  plot([0;X;L],[u0;V;uL],[0;X;L],[u0;Vexact;uL],'--','linewidth',2 )
  axis([0 L  0 2])
 
   pause(0.05)
% determine the maximum relative error
   
   Err(n) =max(abs( (V-Vexact)./Vexact));
   T(n) = t;
end
  figure(2)
  plot(T,Err,'linewidth',2)
  xlabel('t','fontsize',20)
  ylabel('Relative Error','fontsize',20)

