diff --git a/SRI_matlab/main.m b/SRI_matlab/main.m
index 9064c7e25260df116b4dc204f9e49ee860ddc27e..4d7bbf1f2f0832d3c7cabccc3958428ccfcb856f 100644
--- a/SRI_matlab/main.m
+++ b/SRI_matlab/main.m
@@ -1,24 +1,26 @@
 clear figure
 
+N = 180;
+
 S = 990;
-I = 50;
+I = 10;
 R = 0;
-n_pop = S + I + R;
 
-input = [S, I, R];
+% input = [S, I, R];
+input = build_input(S, I, R, 12);
 
 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;
+coeffs.InfectionRate = 0.6;
+coeffs.ImmunisationRate = 0.04;
+coeffs.ImmunisationLossRate = 2.7e-4;
+coeffs.DeathRate = 1.3e-3;
+coeffs.NaturalDeathRate = 2e-5;
+coeffs.NatalityRate = 5e-5;
+% VaccinationRate = build_vacc_rate(0, 2.4e-4, 50, N);
+VaccinationRate = build_vacc_rate(0, 0, 50, N);
 
 Slist = zeros(1, N);
+Wlist = zeros(1, N);
 Ilist = zeros(1, N);
 Rlist = zeros(1, N);
 
@@ -28,7 +30,8 @@ for i=t
     Slist(i) = input(1);
     Ilist(i) = input(2);
     Rlist(i) = input(3);
-    input = step(input, coeffs);
+    Wlist(i) = sum(input(4:length(input)));
+    input = step(input, coeffs, VaccinationRate(i));
 end
 
 % figure
@@ -37,10 +40,31 @@ end
 % plot(t, Ilist);
 % plot(t, Rlist);
 % plot(t, Dlist);
-area([Slist.', Ilist.', Rlist.']);
+area([Slist.', Wlist.', Ilist.', Rlist.']);
 
-legend('Susceptible', 'Infected', 'Immunized');
+legend('Susceptible', 'Waiting', 'Infected', 'Immunized');
 title('Evolution of virus in population');
 xlabel('Time');
 ylabel('Population');
-axis([1, N, 0, inf]);
\ No newline at end of file
+axis([1, N, 0, inf]);
+
+% utility functions
+
+function input=build_input(S, I, R, Waiting)
+    input(1) = S;
+    input(2) = I;
+    input(3) = R;
+    for i = 1:Waiting
+        input(3+i) = 0;
+    end
+end
+
+function vacc_rate=build_vacc_rate(start_rate, end_rate, change_step, N)
+    vacc_rate = zeros(1, N);
+    for i=1:change_step
+        vacc_rate(i) = start_rate;
+    end
+    for i=change_step+1:N
+        vacc_rate(i) = end_rate;
+    end
+end
\ No newline at end of file
diff --git a/SRI_matlab/options.m b/SRI_matlab/options.m
index 3dfd8dc5ffd9fc7c2f1ee3b4fccc6f45a6793595..76d8ad7bb3e235616ec29c6f97b375db8dd9bd50 100644
--- a/SRI_matlab/options.m
+++ b/SRI_matlab/options.m
@@ -3,7 +3,6 @@ classdef options
         InfectionRate
         ImmunisationRate
         ImmunisationLossRate
-        VaccinationRate
         DeathRate
         NaturalDeathRate
         NatalityRate
diff --git a/SRI_matlab/step.m b/SRI_matlab/step.m
index 8e1f987ace602053be54af608cdcf175bf681a0d..55f2f964b8cfbd27cae96501689677a96daf9630 100644
--- a/SRI_matlab/step.m
+++ b/SRI_matlab/step.m
@@ -1,13 +1,28 @@
-function res = step(X, options)
+function res = step(X, options, VaccinationRate)
 S = X(1);
 I = X(2);
 R = X(3);
 
-N = S + I + R;
+N = sum(X);
 
-dS = -options.InfectionRate * S * I/N - options.VaccinationRate * S + options.ImmunisationLossRate * R + options.NatalityRate * N - options.NaturalDeathRate * S;
-dI = options.InfectionRate * S * I/N - options.ImmunisationRate * I - options.DeathRate * I;
-dR = options.ImmunisationRate * I + options.VaccinationRate * S - options.ImmunisationLossRate * R - options.NaturalDeathRate * R;
+new_waiting = zeros(1, length(X) - 3);
+for i=5:length(X)
+    new_waiting(i - 3) = X(i - 1) * (1 - options.NaturalDeathRate);
+end
 
-res = [max(S + dS, 0), max(I + dI, 0), max(R + dR, 0)];
+if length(X) > 3
+    infecter = sum(X(4:length(X)) + I);
+    new_waiting(1) = options.InfectionRate * S * infecter/N;
+    dS = -options.InfectionRate * S * infecter/N;
+    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;
+dR = options.ImmunisationRate * I + VaccinationRate * S - options.ImmunisationLossRate * R - options.NaturalDeathRate * R;
+
+res = [max(S + dS, 0), max(I + dI, 0), max(R + dR, 0), new_waiting];
 end
\ No newline at end of file