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


global B mu beta gamma alpha

B = 0.1;
mu = 0.1;
beta = 0.01;
gamma = 0.1;
alpha = 0.1;


% threshold S = (mu+nu)/beta


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

t0 = 0;
tfinal = 100;


S_ini = [5 10 15 20 25 30];
I_ini = [1 1 1 1 1 1];
V_ini = [1 1 1 1 1 1];

for ind = 1:6

y_ini = [S_ini(ind) I_ini(ind) V_ini(ind)];

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

% plot I versus S

figure(1);hold on; plot(y(:,1),y(:,2));
xlabel('S')
ylabel('I')
hold on; plot(S_ini, I_ini, 'o')



% plot S, I, V versus T

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

figure(3); hold on; plot3(y(:,1),y(:,2),y(:,3));
xlabel('S')
ylabel('I')
zlabel('V')
hold on; plot3(S_ini, I_ini, V_ini, 'o')
box on








end