step.m

function res = step(X, options, VaccinationRate)
S = X(1);
I = X(2);
Q = X(3);
R = X(4);

N = sum(X);

new_waiting = zeros(1, length(X) - 4);
for i=6:length(X)
    new_waiting(i - 4) = X(i - 1) * (1 - options.NaturalDeathRate);
end

if length(X) > 4
    infection = S * (1 - 25 * options.InfectionRate / N) ^ I;
    % infecter = sum(X(5:length(X)) + I);
    % infection = options.InfectionRate * S * infecter/N;
    new_waiting(1) = infection;
    dS = -infection;
    dI = X(length(X));
else
    dS = -options.InfectionRate * S * I/N;
    dI = options.InfectionRate * S * I/N;
end

dS = dS - VaccinationRate * S + options.ImmunisationLossRate * R + options.NatalityRate * N - options.NaturalDeathRate * S;
dI = dI - options.ImmunisationRate * I - options.DeathRate * I - options.QuarantineRate * I;
dQ = -options.ImmunisationRate * Q - options.DeathRate * Q + options.QuarantineRate * I;
dR = options.ImmunisationRate * (I + Q) + VaccinationRate * S - options.ImmunisationLossRate * R - options.NaturalDeathRate * R;

res = [max(S + dS, 0), max(I + dI, 0), max(Q + dQ, 0), max(R + dR, 0), new_waiting];
end