% In this m-file we estimate parameters in an SIR model for fitting
% the polio epidemic data below.
%
%  dS/dt = -a*I*S
%  dI/dt =  a*I*S-b*I
%  dR/dt =  b*I
%  S(0)=S_0, I(0)=I_0, R(0)=R_0.
%
%--------------------------------------------------------------------------
clear all, close all
% Input data and other problem dependent functions and parameters.
yvec=[494 759 1016 1215 1619 2964 8489 22377 32618 38153 41462 42375]';
tvec=[1:1:length(yvec)]';
%        a     b   
b_0 = [0.0002 0.1]'; 
global initial_conditions
initial_conditions = [50831;1;494];
nparams = length(b_0);
npts    = length(yvec);
nlin_fn = @poliofun;

%--------------------------------------------------------------------------
% Compute optimal parameter values using nlinfit and plot.
%--------------------------------------------------------------------------
[b_opt,r,J,C,mse] = nlinfit(tvec,yvec,nlin_fn,b_0);
[tvec,model_opt] = ode45(@polioode,tvec,initial_conditions,[],b_opt);
figure(1)
  subplot(311), plot(tvec,model_opt(:,1))
  subplot(312), plot(tvec,model_opt(:,2))
  subplot(313), plot(tvec,model_opt(:,3),tvec,yvec,'o')