%This is a block of code that solves the model

%These are the initial conditions used in all simulations starting on Feb
%15 2020
%there are 18 compartments
E0 = 33/population;
Ev0 = 0;
Ev20 = 0;
Ev30 = 0;
Ev3O0 = 0;
S0 = (1-E0);
I0 = 0;
Iv0 = 0;
Iv20 = 0;
Iv30 = 0;
Iv3O0 = 0;
R0 = 0;
RO0 = 0;
H0 = 0;
D0 = 0;
CI0 = 0;
V0 = 0;
CItotal0 = 0;


y0 = zeros(17,1);

%there are 24 parameters
params = zeros(25,1);

params(1) = betabar;
params(2) = kappa ;
params(3) = zeta ;
params(4) = gamma;
params(5) = sigma ;
params(6) = eta;
params(7) = etaO;
params(8) = seasonalsize;
params(9) = seasonalposition;
params(10) = fatiguesize;
params(11) = fatiguemean;
params(12) = fatiguesig;
params(13) = betabarv;
params(14) = tv;
params(15) = Evbar;
params(16) = lambda;
params(17) = xi;
params(18) = betabarv2;
params(19) = tv2;
params(20) = Evbar2;
params(21) = betabarv3;
params(22) = tv3;
params(23) = Evbar3;
params(24) = nuO;
params(25) = etaO2;
params(26) = indexparams;
params(27) = xiO;

y0(1) = S0;
y0(2) = E0;
y0(3) = Ev0;
y0(4) = Ev20;
y0(5) = Ev30;
y0(6) = Ev3O0;
y0(7) = I0;
y0(8) = Iv0;
y0(9) = Iv20;
y0(10) = Iv30;
y0(11) = Iv3O0;
y0(12) = R0;
y0(13) = RO0;
y0(14) = H0;
y0(15) = D0;
y0(16) = CI0;
y0(17) = V0;
y0(18) = CItotal0;


opts = odeset('Reltol',1e-13,'AbsTol',1e-14,'Stats','on');
[t1,y1] = ode113(@(t,y) omicronodefileBPEAfastvaccines(t,y,params),[t0 tfinal],y0,opts);


S1 = y1(:,1);
E1 = y1(:,2);
Ev1 = y1(:,3);
Ev21 = y1(:,4);
Ev31 = y1(:,5);
Ev3O1 = y1(:,6);
I1 = y1(:,7);
Iv1 = y1(:,8);
Iv21 = y1(:,9);
Iv31 = y1(:,10);
Iv3O1 = y1(:,11);
R1 = y1(:,12);
RO1 = y1(:,13);
H1 = y1(:,14);
D1 = y1(:,15);
CI1 = y1(:,16);
V1 = y1(:,17);
CItotal1 = y1(:,18);

%this is daily deaths
Ddot1 = zeta*H1;
Ddotdot1 = eta*gamma*(I1+Iv1+Iv21) + etaO*gamma*Iv31 + etaO2*gamma*Iv3O1 - zeta*H1;
gD1 = zeta*Ddotdot1./Ddot1;
