clear all

Population_Total = 15.6*10^6;
Population_Older = 2.78*10^6;
Population_Workers = 0.8*10^6;
Population_Young_NW = Population_Total - Population_Older- Population_Workers;
hospital_rate_O = 125/100000;
hospital_rate_Y = 50/100000;
death_rate_O = 56/100000;
death_rate_Y = 15/100000;
efficacy_O = 56/100;
efficacy_Y = 80/100;
hospital_O(1)=0;
hospital_Y(1)=0;
hospital_W(1)=0;
death_O(1)=0;
death_Y_NW(1)=0;
death_W(1)=0;
death_Y(1) = 0;
dt = 1;
%STRATEGY 1
Old_Sick(1) = 0.3*Population_Older;
Young_Sick_NW(1) = 0.3*Population_Young_NW;
Workers_Sick(1) = 0.3*Population_Workers;
Young_Sick(1) = Young_Sick_NW(1)+Workers_Sick(1)
Old_Hlt(1) = Population_Older - Old_Sick(1);
Young_Hlt_NW(1) = Population_Young_NW - Young_Sick_NW(1);
Workers_Hlt(1)= Population_Workers - Workers_Sick(1);
Young_Hlt(1) = Young_Hlt_NW(1) + Workers_Hlt(1);
k_1 = hospital_rate_O*(1-efficacy_O);
d_1 = death_rate_O*(1-efficacy_O);
k_2 = hospital_rate_Y*(1-efficacy_Y);
d_2 = death_rate_Y*(1-efficacy_Y);
T(1)=0;
t = 1;
while Old_Sick(t) > 0
    hospital_O(t+1) = hospital_O(t)+ k_1*Old_Sick(t)*dt;
    hospital_Y(t+1) = hospital_Y(t)+hospital_rate_Y*Young_Sick_NW(t)*dt;
    hospital_W(t+1) = hospital_W(t)+k_2*Workers_Sick(t)*dt;
    death_O(t+1) = death_O(t) + d_1*Old_Sick(t)*dt;
    death_Y_NW(t+1) = death_Y_NW(t) + death_rate_Y*Young_Sick_NW(t)*dt;
    death_W(t+1) = death_W(t) + d_2*Workers_Sick(t)*dt;
    death_Y(t+1) = death_Y_NW(t+1) + death_W(t+1);
    Old_Sick(t+1)= Old_Sick(t)- death_O(t+1) - hospital_O(t+1);
    Young_Sick_NW(t+1)= Young_Sick_NW(t)- death_Y_NW(t+1) - hospital_Y(t+1);
    Workers_Sick(t+1)= Workers_Sick(t)- death_W(t+1) - hospital_W(t+1);
    Young_Sick(t+1) = Young_Sick_NW(t+1) + Workers_Sick(t+1);
    Old_Hlt(t+1) = Population_Older - Old_Sick(t+1);
    Young_Hlt_NW(t+1) = Population_Young_NW - Young_Sick_NW(t+1);
    Workers_Hlt(t+1)= Population_Workers - Workers_Sick(t+1);
    Young_Hlt(t+1) = Young_Hlt_NW(t+1) + Workers_Hlt(t+1);
    T(t+1) = T(t) +1;
    t=t+1;
end
% plot(T,Old_Sick)
% hold on
% plot(T,Old_Hlt - Old_Hlt(1))
% hold on
% plot(T,death_O)
% hold on
plot(T,Young_Hlt)
hold on
plot(T,death_Y)
hold on
plot(T,Young_Sick)