 % This file solves the diffusion equation using backward Euler method
clear
 
L = 1; %length of the domain
D = 1; % the diffusion coefficient
del0 = 0; % porosity or reaction rate at x=0 (for Robin condition)
delL = 0; % porosity or reaction rate at x=L (for Robin condition)
% Robin BC at x=0:  -Du_x + del0 u0 = U0 at x=0
%Robin BC at x=L: Du_x +delL u_N = UL  
U0 =  0;
UL = 0;

N = 20; %number of grid points = N+1     
 dx = L/(N); % grid size
 
 rsc0 = 2*del0*dx/D; %the robin condition diagonal factor 
 rscL = 2*delL*dx/D;  %the robin condition diagonal factor
 sc =  2*ones(N+1,1) + [rsc0;zeros(N-1,1);rscL]; % diagonal of difference operator
 scu = [2;ones(N-1,1)]; %upper diagonal for difference operator 
 scl = [ones(N-1,1);2]; %lower diagonal for difference operator
  bc = 2*dx/D*[U0;zeros(N-1,1);UL]; % for inhomogeneous BC's
  
M = 20;  % number of time steps
%pick a time step:
dt = dx/5; % time step  %this is much bigger than for forward euler!
scal = D*dt/dx^2;

X = dx*(0:N)';
 
 %V = exp(-(X-dx*N/2).^2./(2*dx));  % initial data is gaussian
  V = 1-cos(pi*X);  % initial data
 Err(1) = 0;
 T(1) = 0;
kstep = 1;
t=0;

% uses backward Euler 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 = Am\(V+scal*bc);
 
t=t + dt;
end
% exact solution
Vexact = 1-exp(-pi^2*t)*cos(pi*X);
figure(1)
  plot(X,V,X,Vexact,'--','linewidth',2 )
  axis([0 L  0 2])
 xlabel('x','fontsize',20)
 ylabel('u(x,t)','fontsize',20)
   pause(0.05)

% determine the maximum relative error
   
   Err(n) =max(abs( (V-Vexact)./Vexact));
   T(n) = t;
end
  figure(2)
  plot(T,Err)

