% matlab code; question marks ? varies according to simulation requirements

for m=1:50 %number of simulations
    
h=1/3; %timestep. rates are in per days, so here timestep h = 8h.
finalday = 250;
numsteps = 1:h:finalday; %create vector for each step, to be used in loop

%preallocate solution vectors for faster computation
S = zeros(length(numsteps),1); %susceptible
E = zeros(length(numsteps),1); %exposed
Ih = zeros(length(numsteps),1); %infected, hospitalised
Ic = zeros(length(numsteps),1); %infected, community
Dh = zeros(length(numsteps),1); %dead, hospitalised
Dc = zeros(length(numsteps),1); %dead, community
Rh = zeros(length(numsteps),1); %recovered, hospitalised
Rc = zeros(length(numsteps),1); %recovered, community

S(1) = ?; %initial conditions
E(1) = 0;
Ih(1) = 0;
Ic(1) = ?;
Dh(1) = 0;
Dc(1) = 0;
Rh(1) = 0;
Rc(1) = 0;

rho = ?; %proportion of serious cases of exposed tested and hospitalised
c = 0.05; %proportion of infected that dies

N = S(1) + E(1) + Ih(1) + Ic(1) + Rh(1) + Rc(1); %effective total population N. changes every loop to exclude D

%set up ODEs
dSdt = @(S,Ic,betaa)- betaa*Ic*S/N;
dEdt = @(E,Ic,S,alpha1,betaa) betaa*Ic*S/N - alpha1*E ;
dIhdt = @(Ih,E,alpha1,gamma1) alpha1*rho*E - gamma1*Ih;
dIcdt = @(Ic,E,alpha1,gamma1) alpha1*(1-rho)*E - gamma1*Ic;
dDhdt = @(Ih,gamma1) gamma1*c*Ih;
dDcdt = @(Ic,gamma1) gamma1*c*Ic;
dRhdt = @(Ih,gamma1) gamma1*(1-c)*Ih;
dRcdt = @(Ic,gamma1) gamma1*(1-c)*Ic;

numalpha = 100; %number of 1/alpha and 1/gamma sampled from Erlang distribution
alpha = zeros(numalpha,1); %note: it's actually 1/alpha and 1/gamma here, we invert later
gamma = zeros(numalpha,1);

%preallocate solution vector for each alpha sample for faster computation
preS = zeros(numalpha,1);
preE = zeros(numalpha,1);
preIh = zeros(numalpha,1);
preIc = zeros(numalpha,1);
preDh = zeros(numalpha,1);
preDc = zeros(numalpha,1);
preRh = zeros(numalpha,1);
preRc = zeros(numalpha,1);

%set values for R for the 3 phases
for k=1:?
    R(k) = ?;
end

for k=?:?
    R(k) = ?;
end

for k=?:(length(numsteps)-1)
    R(k) = ?;
end

%start iterations for each timestep
for i=1:(length(numsteps)-1)
    
    %sample 1/alpha and 1/gamma from Erlang distribution
    for j = 1:numalpha
        alpha(j) = max(1,gamrnd(2,5.2/2));
        gamma(j) = max(1,gamrnd(2,3/2));
    end
    
    %invert to get alpha, gamma
    alpha = 1./alpha; 
    gamma = 1./gamma;    
    
    %4th order Runge Kutta for each alpha and gamma; averaged at the end
    for l = 1:numalpha

        betaa = gamma(l)*R(i);
 
        S1 = h*dSdt(S(i),Ic(i),betaa);    
        E1 = h*dEdt(E(i),Ic(i),S(i),alpha(l),betaa);
        Ih1 = h*dIhdt(Ih(i),E(i),alpha(l),gamma(l));
        Ic1 = h*dIcdt(Ic(i),E(i),alpha(l),gamma(l));
        Dh1 = h*dDhdt(Ih(i),gamma(l));
        Dc1 = h*dDcdt(Ic(i),gamma(l));
        Rh1 = h*dRhdt(Ih(i),gamma(l));
        Rc1 = h*dRcdt(Ic(i),gamma(l));

        S2 = h*dSdt(S(i)+S1/2,Ic(i)+Ic1/2,betaa);      
        E2 = h*dEdt(E(i)+E1/2,Ic(i)+Ic1/2,S(i)+S1/2,alpha(l),betaa);
        Ih2 = h*dIhdt(Ih(i)+Ih1/2,E(i)+E1/2,alpha(l),gamma(l));
        Ic2 = h*dIcdt(Ic(i)+Ic1/2,E(i)+E1/2,alpha(l),gamma(l));
        Dh2 = h*dDhdt(Ih(i)+Ih1/2,gamma(l));
        Dc2 = h*dDcdt(Ic(i)+Ic1/2,gamma(l));
        Rh2 = h*dRhdt(Ih(i)+Ih1/2,gamma(l));
        Rc2 = h*dRcdt(Ic(i)+Ic1/2,gamma(l));
    
        S3 = h*dSdt(S(i)+S2/2,Ic(i)+Ic2/2,betaa);    
        E3 = h*dEdt(E(i)+E2/2,Ic(i)+Ic2/2,S(i)+S2/2,alpha(l),betaa);
        Ih3 = h*dIhdt(Ih(i)+Ih2/2,E(i)+E2/2,alpha(l),gamma(l));
        Ic3 = h*dIcdt(Ic(i)+Ic2/2,E(i)+E2/2,alpha(l),gamma(l));
        Dh3 = h*dDhdt(Ih(i)+Ih2/2,gamma(l));
        Dc3 = h*dDcdt(Ic(i)+Ic2/2,gamma(l));
        Rh3 = h*dRhdt(Ih(i)+Ih2/2,gamma(l)); 
        Rc3 = h*dRcdt(Ic(i)+Ic2/2,gamma(l));
    
        S4 = h*dSdt(S(i)+S3,Ic(i)+Ic3,betaa);    
        E4 = h*dEdt(E(i)+E3,Ic(i)+Ic3,S(i)+S3,alpha(l),betaa);
        Ih4 = h*dIhdt(Ih(i)+Ih3,E(i)+E3,alpha(l),gamma(l));
        Ic4 = h*dIcdt(Ic(i)+Ic3,E(i)+E3,alpha(l),gamma(l));
        Dh4 = h*dDhdt(Ih(i)+Ih3,gamma(l));
        Dc4 = h*dDcdt(Ic(i)+Ic3,gamma(l));
        Rh4 = h*dRhdt(Ih(i)+Ih3,gamma(l));
        Rc4 = h*dRcdt(Ic(i)+Ic3,gamma(l));
    
        preS(l) = S(i) + (1/6)*(S1+2*S2+2*S3+S4);
        preE(l) = E(i) + (1/6)*(E1+2*E2+2*E3+E4);
        preIh(l) = Ih(i) + (1/6)*(Ih1+2*Ih2+2*Ih3+Ih4);
        preIc(l) = Ic(i) + (1/6)*(Ic1+2*Ic2+2*Ic3+Ic4);
        preDh(l) = Dh(i) + (1/6)*(Dh1+2*Dh2+2*Dh3+Dh4);
        preDc(l) = Dc(i) + (1/6)*(Dc1+2*Dc2+2*Dc3+Dc4);
        preRh(l) = Rh(i) + (1/6)*(Rh1+2*Rh2+2*Rh3+Rh4);
        preRc(l) = Rc(i) + (1/6)*(Rc1+2*Rc2+2*Rc3+Rc4);
    
    end
    
    %average
    S(i+1) = mean(preS);
    E(i+1) = mean(preE);
    Ih(i+1) = mean(preIh);
    Ic(i+1) = mean(preIc);
    Dh(i+1) = mean(preDh);
    Dc(i+1) = mean(preDc);
    Rh(i+1) = mean(preRh);
    Rc(i+1) = mean(preRc);
    
    %update N
    N = S(i+1) + E(i+1) + Ih(i+1) + Ic(i+1) + Rh(i+1) + Rc(i+1);
end

%plot our simulation
p=plot(Ic+Ih,'Color', [0.0 0.0 0.0])
xticks([0 90 180 270 360 450 540 630 720])

%optional here: plot data points for comparison

end

%label axes
xlabel('Time (8 Hours Per Step)')
ylabel('?')