clear figure

S = 990;
I = 50;
R = 0;
n_pop = S + I + R;

input = [S, I, R];

coeffs = options;
coeffs.InfectionRate = 0.09;
coeffs.ImmunisationRate = 0.005;
coeffs.ImmunisationLossRate = 0.001;
coeffs.VaccinationRate = 0.02;
coeffs.DeathRate = 0.04;
coeffs.NaturalDeathRate = 0.002;
coeffs.NatalityRate = 0.003;

N = 1000;

Slist = zeros(1, N);
Ilist = zeros(1, N);
Rlist = zeros(1, N);

t = 1:N;

for i=t
    Slist(i) = input(1);
    Ilist(i) = input(2);
    Rlist(i) = input(3);
    input = step(input, coeffs);
end

% figure
% hold on
% plot(t, Slist);
% plot(t, Ilist);
% plot(t, Rlist);
% plot(t, Dlist);
area([Slist.', Ilist.', Rlist.']);

legend('Susceptible', 'Infected', 'Immunized');
title('Evolution of virus in population');
xlabel('Time');
ylabel('Population');
axis([1, N, 0, inf]);