clear all

  global N
 
  S0 = 900;
  I0 = 100;
  R0 = 0;
  N = 1000;

  tf = 20;
  dt = 0.05;
  t_data = 0:dt:tf;
  ind = length(t_data);
  Y0 = [S0; I0; R0];
  I_mu = 0.5;
  I_eta = 0.1;
  I_gamma = 1.5;

  alpha = 0.2;
  beta = 15;
%  alpha = 5;
%  beta = 9;

  ode_options = odeset('RelTol',1e-6);

  r = 1000;
  ell = 1000;
  Delta = 1/ell;

  tic
  wait = waitbar(0,'Please wait...');
  for i = 1:r
    %params_star = rand(1,4);
    params_star(1) = I_eta*rand(1,1); 
    params_star(2) = betarnd(alpha,beta,1);
    params_star(3) = I_gamma*rand(1,1);
    params_star(4) = I_mu*rand(1,1);
  
    for j=1:4
      params = params_star;
      params(j) = params_star(j) + Delta; 
      [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
      y_A(j) = sum(dt*Y(:,3));
    end
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params_star);
    y_star = sum(dt*Y(:,3));
 
    d1(i) = (y_A(1) - y_star)/Delta;
    d2(i) = (y_A(2) - y_star)/Delta;
    d3(i) = (y_A(3) - y_star)/Delta;
    d4(i) = (y_A(4) - y_star)/Delta;
    waitbar(i/r);
  end

  mu1 = (1/r)*sum(d1);
  mu2 = (1/r)*sum(d2);
  mu3 = (1/r)*sum(d3);
  mu4 = (1/r)*sum(d4);

  mu1s = (1/r)*sum(abs(d1));
  mu2s = (1/r)*sum(abs(d2));
  mu3s = (1/r)*sum(abs(d3));
  mu4s = (1/r)*sum(abs(d4));
  mus = [mu1s mu2s mu3s mu4s]

  sigma1 = (1/(r-1))*sum((d1-mu1).^2);
  sigma2 = (1/(r-1))*sum((d2-mu2).^2);
  sigma3 = (1/(r-1))*sum((d3-mu3).^2);
  sigma4 = (1/(r-1))*sum((d4-mu4).^2);
  sigma2 = [sigma1 sigma2 sigma3 sigma4]
  sigma = sqrt(sigma2)