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.2;

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

  A(:,1) = I_eta*rand(M,1);
  A(:,2) = betarnd(alpha,beta,M,1);
  A(:,3) = I_gamma*rand(M,1);
  A(:,4) = I_mu*rand(M,1);
  B(:,1) = I_eta*rand(M,1);
  B(:,2) = betarnd(alpha,beta,M,1);
  B(:,3) = I_gamma*rand(M,1);
  B(:,4) = I_mu*rand(M,1);
  C1 = A;
  C2 = A;
  C3 = A;
  C4 = A;
  C1(:,1) = B(:,1);
  C2(:,2) = B(:,2);
  C3(:,3) = B(:,3);
  C4(:,4) = B(:,4);

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

  tic
  wait = waitbar(0,'Please wait...');
  for j=1:M
    params = A(j,:);
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
    Y_A(j,:) = Y(:,3)';

    params = B(j,:);
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
    Y_B(j,:) = Y(:,3)';

    params = C1(j,:);
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
    Y_C1(j,:) = Y(:,3)';

    params = C2(j,:);
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
    Y_C2(j,:) = Y(:,3)';

    params = C3(j,:);
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
    Y_C3(j,:) = Y(:,3)';

    params = C4(j,:);
    [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
    Y_C4(j,:) = Y(:,3)';

    waitbar(j/M)
  end
  toc

  Y_D = [Y_A; Y_B];

  num1 = 0;
  num2 = 0;
  num3 = 0;
  num4 = 0;
  Num1 = 0;
  Num2 = 0;
  Num3 = 0;
  Num4 = 0;
  den = 0;
  for i = 1:ind
    y_A = Y_A(:,i);
    y_B = Y_B(:,i);
    y_C1 = Y_C1(:,i);
    y_C2 = Y_C2(:,i);
    y_C3 = Y_C3(:,i);
    y_C4 = Y_C4(:,i);
    y_D = Y_D(:,i);

    f02 = ((1/(2*M))^2)*sum(y_D)*sum(y_D);
    den = den + dt*((1/(2*M))*y_D'*y_D - f02);

    num1 = num1 + dt*((1/M)*(y_B'*y_C1 - y_B'*y_A));
    num2 = num2 + dt*((1/M)*(y_B'*y_C2 - y_B'*y_A));
    num3 = num3 + dt*((1/M)*(y_B'*y_C3 - y_B'*y_A));
    num4 = num4 + dt*((1/M)*(y_B'*y_C4 - y_B'*y_A));

    Num1 = Num1 + dt*((1/(2*M))*(y_A'*y_A - 2*y_A'*y_C1 + y_C1'*y_C1));
    Num2 = Num2 + dt*((1/(2*M))*(y_A'*y_A - 2*y_A'*y_C2 + y_C2'*y_C2));
    Num3 = Num3 + dt*((1/(2*M))*(y_A'*y_A - 2*y_A'*y_C3 + y_C3'*y_C3));
    Num4 = Num4 + dt*((1/(2*M))*(y_A'*y_A - 2*y_A'*y_C4 + y_C4'*y_C4));

    S1(i,1) = num1/den;
    S2(i,1) = num2/den;
    S3(i,1) = num3/den;
    S4(i,1) = num4/den;

    ST1(i,1) = Num1/den;
    ST2(i,1) = Num2/den;
    ST3(i,1) = Num3/den;
    ST4(i,1) = Num4/den;

  end

  S = [S1 S2 S3 S4];
  ST = [ST1 ST2 ST3 ST4];

%
%  Plot the results
%

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

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

  figure(2)
  plot(t,S(:,1),'k',t,S(:,2),'b--',t,S(:,3),'r-.',t,S(:,4),'k:','linewidth',3)
  axis([0 20 0 0.9])
  set(gca,'Fontsize',[22]);
  legend('\eta','k','\gamma','\mu','Location','NorthEast')
  xlabel('Time (Days)')
  ylabel('Generalized Sobol'' Indices')
  %print -depsc Gen_Sobol_Beta2_15
  %print -depsc Gen_Sobol_Beta5_9

  figure(3)
  plot(t,ST(:,1),'k',t,ST(:,2),'b--',t,ST(:,3),'r-.',t,ST(:,4),'k:','linewidth',3)
  axis([0 20 0 0.9])
  set(gca,'Fontsize',[22]);
  legend('\eta','k','\gamma','\mu','Location','NorthEast')
  xlabel('Time (Days)')
  ylabel('Generalized Total Sobol'' Indices')
  %print -depsc Gen_Total_Sobol_Beta2_15
  %print -depsc Gen_Total_Sobol_Beta5_9