xlab = 'Week ($n$)'; 
ylab = 'Infected Cases ($I_n$)'; 		
% mytitle = ''; 	% Title
load('sir_data');
tt = sirn;
popi = siri;
K = length(sirn);
N = 157759;
p = [3.992759714, 3.517005231];
Sn(1) = [157756];
In(1) = [3];
Rn(1) = [0];
for i = 2:K
    Sn(i) = Sn(i-1) - p(1)*Sn(i-1)*In(i-1)/N;
    In(i) = In(i-1) + p(1)*Sn(i-1)*In(i-1)/N - p(2)*In(i-1);
    Rn(i) = Rn(i-1) + p(2)*In(i-1);
end
IVn(1) = [3];
SVn(1) = [0.95*N - IVn(1)];
RVn(1) = [0.05*N];
for i = 2:K
    SVn(i) = SVn(i-1) - p(1)*SVn(i-1)*IVn(i-1)/N;
    IVn(i) = IVn(i-1) + p(1)*SVn(i-1)*IVn(i-1)/N - p(2)*IVn(i-1);
    RVn(i) = RVn(i-1) + p(2)*IVn(i-1);
end
IQn(1) = [3];
SQn(1) = [N - IQn(1)];
RQn(1) = [0];
pq(1) = 0.95*p(1);
for i = 2:K
    SQn(i) = SQn(i-1) - pq(1)*SQn(i-1)*IQn(i-1)/N;
    IQn(i) = IQn(i-1) + pq(1)*SQn(i-1)*IQn(i-1)/N - p(2)*IQn(i-1);
    RQn(i) = RQn(i-1) + p(2)*IQn(i-1);
end
ITn(1) = [3];
STn(1) = [N - ITn(1)];
RTn(1) = [0];
pt(2) = 1.05*p(2);
for i = 2:K
    STn(i) = STn(i-1) - p(1)*STn(i-1)*ITn(i-1)/N;
    ITn(i) = ITn(i-1) + p(1)*STn(i-1)*ITn(i-1)/N - pt(2)*ITn(i-1);
    RTn(i) = RTn(i-1) + pt(2)*ITn(i-1);
end
plot(tt,popi,'bo','MarkerSize',8,'LineWidth',1.5);
hold on
plot(tt,In,'r-','LineWidth',1.5);
plot(tt,IVn,'-','color',[0,0.5,0],'LineWidth',1.5);
plot(tt,IQn,'m-','LineWidth',1.5);
plot(tt,ITn,'k-','LineWidth',1.5);

grid
legend('CDC Data','SIR Model','Vaccine','Quarantine',...
    'Tamiflu','Location','Northwest');
xlim([0 48]);
ylim([0 1500]);
fontlabs = 'Times New Roman';   % Font type used in labels
xlabel(xlab,'FontSize',14,'FontName',fontlabs,'interpreter','latex');  
ylabel(ylab,'FontSize',14,'FontName',fontlabs,'interpreter','latex'); 
%title(mytitle,'FontSize',16,'FontName',fontlabs,'interpreter','latex'); 
set(gca,'FontSize',12); 
print -depsc sir_treat.eps
