%% bet= 2.0  -- 1/mu = 5 days
%% integrate for 75 days.
%% Parameter R0 = bet* N / mu -- N = N susceptible hosts
%% I = infectious --> transmits virus
%% S = Susceptible --> not immune and susceptible of getting infected
%% R = 'removed' - either recovered or dead.
%% e.g.  N = 100,000  I = 1000 R =0 as initial state. 

%% reproduction number =
%% do it for R = 3. plot, then R = 1.5 replot and then R = 1.1 replot

%%
%%====================
%% total pop. N = 100,000 - Initially: S = 99,000 . I = 1000 and R = 0;
N = 100000;
I = 1000;
S = N - I;
R = 0;
%%-------------------- transmission rates etc
%% R_0 = 3.0 then 1.7 -- then 1.
mu = 1/5;    %% infection period 
R0 = 1.5;
bet = mu * R0/N

%%====================integration scheme:
%%                    forwared euler

Y = [S; I; R]
Dt = 0.5;
It = Y;
days(1) = 0.0;

for j=2:150
  days(j) = days(j-1)+ Dt;
  M = [-bet*Y(2), 0, 0;
       0, bet*Y(1)-mu, 0;
       0, mu, 0];
  Y = Y + Dt*M*Y;
  It(:,j) = Y;
end
%%plot(days,0.05*It(3,:),'k-','linewidth',2)
plot(days,It(2,:),'r--','linewidth',2)

hold on

%%h=legend('R0=3.0','R0=1.4','R0=1.1','R0=1.0');
%%set(h,'fontsize',18)

xlabel('Days')
ylabel('Number infected (I)')
