 % This file solves the diffusion equation using Crank Nicholson
 % can be used for Neumann and Robin conditions
clear
 
% this code is used to solve one of the exercises 
L = 1; %length of the domain
D = 1; % the diffusion coefficient
del0 = 0; % porosity or reaction rate at x=0 (for Robin condition)
delL = 1; % 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 actual grid points = N+1 
% N-1 interior grid points
% there are N+1 unknowns
 dx = L/N; % grid size
 
 rsc0 = 2*del0*dx/D; %the robin condition diagonal factor at x=0
 rscL = 2*delL*dx/D;  %the robin condition diagonal factor at x=L
 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
    
%pick a time step:
  
M = 100;  % 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

X = dx*(0:N)';
 
%specify initial data
  
 V = ones(N+1,1);
 Err(1) = 0;
 T(1) = 0;
 V0(1) = 0;
 
kstep = 5;
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*bc);
V = Am\(V+V_t);
 
t=t + dt;
end
% exact solution
 T(n) = t;
 V0(n) = V(1);
 
figure(1)
  plot(X,V,'linewidth',2 )
  axis([0 L  0 2])
 xlabel('x','fontsize',20)
 ylabel('u(x,t)','fontsize',20)
   pause(0.1)

   
   % determine the maximum relative error
   
   
end

p = polyfit(T(10:end),log(V0(10:end)),1)

  figure(2)
  plot(T,log(V0),T,p(2)+p(1)*T,'--','linewidth',2)

  % estimate the root of the equation a=1=m*tan(m);
  m = [0:.001:2];
  a = m.*tan(m);
  mdx =  min(find(a>=1));
  ms = m(mdx-1)+(1-a(mdx-1))*(m(mdx)-m(mdx-1))/(a(mdx)-a(mdx-1));
  
  lambda = ms^2
  axis([0 2 0 2])