 % This file solves the diffusion equation using Crank Nicholson
 % can be used for Neumann and Robin conditions
clear
 
 
L = 1; %length of the domain
D = 1; % the diffusion coefficient
del0 = 1; % 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 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
    
%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

X = dx*(0:N)';
 
%specify initial data
V=ones(N+1,1);  % uniform initial data

 
kstep = 5;
t=0;
T(1) = 0;
Vmid(1) = 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
T(n) = t;
Vmid(n) = 1-V(N/2+1);
figure(1)
  plot(X,V,X,1-V,'linewidth',2 )
  axis([0 L  0 1])
  xlabel('\xi','fontsize',20)
 ylabel('u(\xi,\tau), v(\xi,\tau)','fontsize',20)
   pause(0.1)

end

figure(2)
plot(T,Vmid,'linewidth',2)
xlabel('\tau','fontsize',20)
ylabel('V(1/2,\tau)','fontsize' ,20)
  
