function Simulate(DiffStep,dt)
% ----------------------------------------------------------------------------------- 
% Animate the solution to the Advection-Diffusion-Reaction equation using either 
% the TR-BDF2 rule for the Diffusion step (DiffStep = 2, added stability for larger dt)
% or the Crank-Nicholson rule (DiffStep = 1, less stable for larger dt)
% @ Sohan Dharmaraja MIT JUN2007

clc
p=20; m=p^2;
dx = 1/(p+1);
%---------------------------------------------------------------------------------------------------------------

% compute the initial solution profile x0 (soln to heat equation)
% Labelling of the nodes for heat sources is: 
% f(1,1) is the bottom RH element
% f(2,1) is the one to the left of the bottom RH element
% f(p,1) is the bottom LH element
% f(p+1,1) is the element ABOVE f(1,1), etc

Ab=zeros(m,p+1); Ab(:,1)=4; Ab(:,2)=-1; Ab(p:p:m,2)=0; Ab(:,p+1)=-1;
A=spdiags(Ab(:,[1,2,p+1]),-[0,1,p],m,m); A=A+tril(A,-1)';


f=10*ones(m,1)/(p+1)^2;
% f(p*floor(p/2)+ floor(p/2),1)=10/((p+1)^2); % almost centred heat source
%f(1,1)=10/((p+1)^2); f(p,1) = f(1,1); f(m,1) = f(1,1); f(m-p+1,1) = f(1,1); %heat sources at the corners
x0=A\f;

ninit  = x0;
cstring = (dx^3)*(1:p)';
ctmp = [];
for y = 1:p
    ctmp = cat(1,ctmp,cstring);
end
cinit = ctmp;

nsteps = 1;
Dn     = 0.05;
Cn     = 0.01;

FirstPlot(p,x0)
v = caxis;
close;


ncur = ninit;
ccur = cinit;

% Plot the initial n profile
% ------------------------------------------------------------------
PlotData(p,ncur,v)
pause(2);
% ------------------------------------------------------------------

for i=1:1000
        
    ccur = Reaction(cinit,ninit,dt/2,nsteps);
    ncur = Advection(ninit,ccur,p,dt/2,nsteps);

if DiffStep == 1
    ncur = CNDiffusion(Dn,ncur,p,dt,nsteps);
    ccur = CNDiffusion(Cn,ccur,p,dt,nsteps);
else
    ncur = TRBDFDiffusion(Dn,ncur,p,dt,nsteps);
    ccur = TRBDFDiffusion(Cn,ccur,p,dt,nsteps);
end
    ncur = Advection(ncur,ccur,p,dt/2,nsteps);
    ccur = Reaction(ccur,ncur,dt/2,nsteps);
        
if mod(i,10) == 0
        PlotData(p,ncur,v)
        pause(0.3)
end

    cinit = ccur;
    ninit = ncur;
    close;
end