%
%                       SIR_OAT_1D
%

  clear all
  close all

  global N

%
% Set the initial conditions and number in the population.
%
  S0 = 900;
  I0 = 100;
  R0 = 0;
  Y0 = [S0; I0; R0];
  N = 1000;

%
% Construct the time interval.
%

  tf = 5;
  dt = 0.01;
  t_data = 0:dt:tf;
  ind = length(t_data);

%
%  Randomly sample values of gamma, k, r and delta and construct parameter matrices
%

  M = 100;                          % Number of random samples -- you can increase this in your final run

  gamma_nom = 0.9;                  % Nominal value of gamma
  k_nom = 0.9;
  r_nom = 0.9;
  delta_nom = 0.9;

  per = 0.2;                        % Percentage used when computing intervals

  xl = gamma_nom - per*gamma_nom;   % Lower endpiont of interval
  xr = gamma_nom + per*gamma_nom;   % Upper endpoint of interval

  Gamma_vals = (gamma_nom - per*gamma_nom) + 2*per*gamma_nom*rand(M,1);
  K_vals = (k_nom - per*k_nom) + 2*per*k_nom*rand(M,1);
  R_vals = (r_nom - per*r_nom) + 2*per*r_nom*rand(M,1);
  Delta_vals = (delta_nom - per*delta_nom) + 2*per*delta_nom*rand(M,1);

%
%  Params1 contains the values of gamma sampled out of the interval [xl, xr] with fixed
%  values of k, r and delta
%
%  Params contains randomly sampled values of gamma, k, r and delta
%

  Params1 = [Gamma_vals r_nom*ones(M,1) delta_nom*ones(M,1) k_nom*ones(M,1)];
  Params = [Gamma_vals R_vals Delta_vals K_vals];
  
%
% Set the ode options and solve the system M times
%

  ode_options = odeset('RelTol',1e-6);     % Set the ode relative tolerance

  h = waitbar(0,'Please wait...');
  for j=1:M
    params1 = Params1(j,:);
    [t,Y] = ode15s(@SIR_rhs,t_data,Y0,ode_options,params1);
    Response1(j) = sum(dt*Y(:,3));

    params = Params(j,:);
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
    Response(j) = sum(dt*Y(:,3));
    waitbar(j/M)
  end
  close(h)

  Sval = Y(:,1);
  Ival = Y(:,2);
  Rval = Y(:,3);

%
% Plot the solutions.
%

  yl = min(Response);
  yr = max(Response);

  figure(1) 
  plot(t,Sval,'b',t,Ival,'k--',t,Rval,'r-.','linewidth',3)
  set(gca,'Fontsize',[20]);
  legend('Susceptible','Infected','Recovered','Location','Northeast')
  xlabel('Time')
  ylabel('Number of Individuals')


  figure(2)
  subplot(2,2,1)
  plot(Gamma_vals,Response1,'x')
  axis([xl xr yl yr])
  set(gca,'Fontsize',[20]);
  xlabel('\gamma')
  ylabel('y')