
  clear all

  global N

%
% Set the level of interactions and constants in the SIR model.  Also construct the time interval.
%

  alpha = input('alpha = ');
  beta = input('beta = ');
 
  S0 = 900;
  I0 = 100;
  R0 = 0;
  N = 1000;

  tf = 5;
  dt = 0.01;
  t_data = 0:dt:tf;
  ind = length(t_data);
  Y0 = [S0; I0; R0];

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

%
% Specify the number of random evaluation points r and stepsize ell.
%

  r = 20;
  ell = 40;

  val1 = 1;
  val2 = 1;
  val3 = 1;
  val4 = 1;

%
% Compute the Morris indices.
%

  for i = 1:r
    Num = rand(1,4);
    Num(2) = betarnd(alpha,beta,1);
    for j = 1:ell-2
      v1 = abs(j/ell - Num(1));
      v2 = abs(j/ell - Num(2));
      v3 = abs(j/ell - Num(3));
      v4 = abs(j/ell - Num(4));
      if v1 < val1
         q1 = j/ell;
         val1 = v1;
      end
      if v2 < val2
         q2 = j/ell;
         val2 = v2;
      end
      if v3 < val3
         q3 = j/ell;
         val3 = v3;
      end
      if v4 < val4
         q4 = j/ell;
         val4 = v4;
      end
    end
    q = [q1 q2 q3 q4];
    Q(i,:) = q;
%    Delta = ell/(2*(ell-1));
    Delta = 1/ell;
    
    num = rand(1);
    if 0 < num & num < .25
      P = [1 0 0 0;0 1 0 0; 0 0 1 0; 0 0 0 1];
    elseif 0.25 < num & num < 0.5
      P = [0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1];
    elseif 0.5 < num & num < 0.75
      P = [0 0 1 0; 1 0 0 0; 0 1 0 0; 0 0 0 1];
    else
      P = [0 0 0 1; 1 0 0 0; 0 1 0 0; 0 0 1 0];
    end

    for j=1:4
      num = rand(1);
      if num < 0.5
        D(j,j) = -1;
      else
        D(j,j) = 1;
      end
    end
    J1 = ones(5,1);
    J3 = ones(5,4);
    B = [0 0 0 0;1 0 0 0;1 1 0 0; 1 1 1 0; 1 1 1 1];
    Bs = (J1*q + (Delta/2)*((2*B - J3)*D + J3))*P;
  
    for j=1:5
      params = Bs(j,:);
      [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
      y_A(j) = sum(dt*Y(:,3));
    end
 
    for j=1:4
      if Bs(j+1,1)~=Bs(j,1)
        d1(i) = (y_A(j+1) - y_A(j))/Delta;
      elseif Bs(j+1,2)~=Bs(j,2)
        d2(i) = (y_A(j+1) - y_A(j))/Delta;
      elseif Bs(j+1,3)~=Bs(j,3)
        d3(i) = (y_A(j+1) - y_A(j))/Delta;
      elseif Bs(j+1,4)~=Bs(j,4)
        d4(i) = (y_A(j+1) - y_A(j))/Delta;
      end
    end
  end

%
% Compute the Morris indices.
%

  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);
  sigma = [sqrt(sigma1) sqrt(sigma2) sqrt(sigma3) sqrt(sigma4)]
