% This is the main function to solve 
% the SIR model 


global B Bp mu beta mup mubar


B = 0.1;
Bp = 0.1;
mu = 0.1;
beta = 0.5;
mup = 0.1;
mubar = 0.2;

% threshold 
Sbar = B/mu;
Spbar = Bp/mup;

thres = mubar*mup-beta^2*Sbar*Spbar


options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1.e-4 1.e-4]);

t0 = 0;
tfinal = 100;


S_ini = [1 2 3 4 5 6];
I_ini = [1 1 1 1 1 1];
Sp_ini = [1 1 1 1 1 1];
Ip_ini = [1 1 1 1 1 1];

for ind = 1:6

y_ini = [S_ini(ind) I_ini(ind) Sp_ini(ind) Ip_ini(ind)];

[T,y] = ode45(@hum_mos_eqn,[t0 tfinal],y_ini,options);





% plot S, I, Sp, Ip versus T

figure(2);hold on; plot(T,y(:,1),T,y(:,2),T,y(:,3),T,y(:,4));
xlabel('t')
legend('S','I','Sp','Ip')





end