 % This file solves the diffusion equation using Crank Nicholson
 % can be used for Neumann and Robin conditions
clear
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);
 
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 = 200; %number of grid points = N+1 (N needs to be an even number)    
 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 = 15;  % 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)';
cset=[0,0.5];
 for js = 1:2
     cosamp=cset(js)
%specify initial data
 V = 1.5*exp(-(X-dx*N/2).^2./(2*dx))+cosamp*cos(20*pi*X);  % initial data is gaussian
 
 V = 1.5*cos(pi*(X-dx*N/2))+cosamp*cos(20*pi*X);
 figure(js)
  plot(X,V,'--','linewidth',2 )
  axis([0 L  -1 2])
  xlabel('x','fontsize',20)
 ylabel('u(x,t)','fontsize',20)
   pause(0.1)
hold on
Vs = V;
kstep = 2;
t=0;
Ts = t;


% 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;
Ts=[Ts,t];

Vs = [Vs,V];
end
figure(js)
  plot(X,V,'linewidth',2 )
  axis([0 L  -1 2])
  xlabel('x','fontsize',20)
 ylabel('u(x,t)','fontsize',20)
   pause(0.1)
hold on
end
hold off
if (js==1) 
bottom=min(min(Vs));
top=max(max(Vs));
end
[tt,xx]=meshgrid(Ts,X);
figure(2+js)
 
mesh(tt,xx,Vs)
caxis manual
caxis([bottom top]);
set(gca,'Ytick',0:0.25:1)
set(gca,'Xtick',0:0.01:0.04)
xlabel('t','fontsize',20)
 ylabel('x','fontsize',20)
 zlabel('u','fontsize',20)
 axis([0 0.04 0 1 -1 2])
 
 end
  
