% This file is solving the SIR equations with diffusive spread

clear
time=cputime;

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

 
%moviename = ['MOVIE.qt']
%MakeQTMovie('start', moviename)

L = 80; % length of the domain
d = 1; % diffusion coefficient  

N=250;                   % number of spatial grid points
M = 25;  %number of time steps
h = L/N; % spatial discretization
dt = h/2;% time step size
K = 1; %reaction rate = 1 wlog
mu = 0.5; %recovery rate the only (dimensionless) parameter
 
scal = d*dt/h^2;

X = h*(1:N)';
 U = zeros(N,1);  
 
 U =0.1* sech( 2*(X-L/2)).^2;  %initial infectives
 V = 1-U; % susceptibles
  % V+U = 1;
 figure(1)
  subplot(2,1,1)
  plot(X, U,'--','linewidth',2)
  xlabel('x','fontsize',20)
   ylabel('I(x,t)','fontsize',20)
    hold on 
  subplot(2,1,2)
  plot(X, V,'--','linewidth',2)
   xlabel('x','fontsize',20)
    ylabel('S(x,t)','fontsize',20)
  axis([0 L  0 1])
  hold on
 % axis([0 L -.5 1])

kstep = 8;
t=0;

% uses Crank Nicolson to solve the diffusion equation

% set up matrix

Atm=(1+scal)*ones(N,1);
Atm(1,1) = 1+scal/2;
Atm(N,1) = 1+scal/2;

Am=diag(Atm)-diag(scal/2*ones(N-1,1),1)-diag(scal/2*ones(N-1,1),-1);

% The algorithm
 

for n=2:M
for k = 1:kstep
F = K*U.*V;   
G = F -mu*U; % infective reaction rate
U_t =  scal/2 * ([U(2:N );U(N )] - 2*U + [U(1 );U(1:N-1 )]) +dt*G ;

U = Am\(U+U_t);  %update the infectives
V = V -dt*F; % update susceptibles
t=t + dt;
end
figure(1)
  subplot(2,1,1)
  plot(X,U)
  xlabel('\xi','fontsize',20)
   ylabel('u(\xi,\tau)','fontsize',20)
  hold on
  subplot(2,1,2)
  plot(X, V )
   xlabel('\xi','fontsize',20)
  ylabel('\sigma(\xi,\tau)','fontsize',20)
  axis([0 L  0 1])
 % axis([0 L -.5 1])

hold on
%MakeQTMovie('addframe','moviename')
%MakeQTMovie('quality', 1.0)

  pause(0.05)

end
 

%MakeQTMovie('finish','moviename')
  

TL = t;
time=cputime-time        % outputs the cpu time taken to solve the equations

time_steps = M;

