% https://www.nber.org/papers/w26867.pdf
global beta sigma gamma f N

t0=0;
tf=150;
N=3.3E8;
NI=33;
NE=4*NI;
f=0;

Recovery = 11.5;
Incubation  = 5.6;
gamma = 1.0/Recovery;
sigma = 1/Incubation;
R0 = 4.7;
beta = R0*gamma;

[t,y]=ode15s('seirf',[t0,tf],[N-NI-NE,NE,NI,0]);
plot(t,y(:,2),t,y(:,3),t,y(:,4))
xlabel('Time')
ylabel('Populations')
title('SEIR Model')
legend('Exposed','Infectious','Removed','Location','northwest')

% Get Data
M=79;
M0=15;
Brange=['B',num2str(M0),':B',num2str(M)];
Crange=['C',num2str(M0),':C',num2str(M)];
A = xlsread('COVID19data.xlsx', 1, Brange);
C = xlsread('COVID19data.xlsx', 1, Crange);
hold on
%plot(0:M-M0,A,'ro',0:M-M0,C,'bo')
plot(0:M-M0,C,'bo')
hold off
%axis([0 60 0 2e06])