% Driver file that specifies parameter values, calls ode solver and plots semilog graph of results

tmax=200; % length of time (yrs) to integrate
tspan=[0 tmax];
N=1000000; % population size
Sinit=0.1*N; % initial susceptible population
Iinit=1; % initial infected population
Rinit=N-Sinit-Iinit; % initial recovered/immune population

R0=17;
gamma=365/10; % inverse of average infectious period (yrs)
mu=1/70; % inverse of life expectancy (yrs)
beta0=R0*(gamma+mu);
beta1=0.1;

options = odeset('RelTol',1e-6,'AbsTol',1e-12); % solver options
[T,Y] = ode45(@SIRforced,tspan,[Sinit Iinit Rinit],options,beta0,beta1,gamma,mu);% solve odes
plot(T,Y(:,2))
xlim([180 200])
xlabel('Time (years)')
ylabel('I(t)')




