 % This file solves the diffusion equation with Neumann BC's at x=0,
 % Dirichlet BC's at x=L
 %using Crank Nicholson method for an exercise
clear

L = 1; %length of the domain
D = 1; % the diffusion coefficient
del0 = 0; % porosity or reaction rate at x=0 (for Robin condition)

% Robin BC at x=0:  -Du_x + del0 u0 = U0 at x=0

U0 = 0;

% Dirichlet condition at x=L u = UL
UL = 0;

N = 20; 
%number of grid points = N+1, N-1 interior points     
 dx = L/N; % grid size
 % there are N unknowns
 
 rsc0 = 2*del0*dx/D; %the robin condition diagonal factor at x=0
 
 sc =  2*ones(N,1) + [rsc0;zeros(N-1,1)]; % diagonal of difference operator
 scu = [2;ones(N-2,1) ]; %upper diagonal for difference operator 
 scl = [ones(N-1,1)]; %lower diagonal for difference operator
 bc = [2*dx/D*U0;zeros(N-2,1);UL]; % for inhomogeneous BC's
   
 X = [0:N-1]'*dx;  % x values for unknowns

%pick a time step:
  
M = 100;  % number of time steps
 
dt = dx/10; % 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 = ones(N,1);  % initial data for exercise
kstep = 5;
t=0;
T(1) = 0;
V0(1) = 1;
% 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 );0] - sc.*V + [0;scl].*[0;V(1:N-1 )]+2*bc);
 
V = Am\(V+V_t);
 
t=t + dt;
end
T(n) = t;
V0(n) = V(1);
% exact solution
figure(1)
  plot([ X;L],[ V;UL],'linewidth',2 )
  axis([0 L  0 1])
 
   pause(0.05)
end
   
 p=polyfit(T(10:end),log(V0(10:end)),1)
 pi^2/4
figure(2)
plot(T,log(V0),T,p(2) + p(1)*T,'--','linewidth',2)

xlabel('t','fontsize',20)
ylabel('ln(V(0,t))','fontsize',20)
