function [G,I,B,R,T]=fingIPVint(tf)

% This function uses matlab's ode45 to numerically solve the IPV system in
% the paper by Eichner and Hadeler (see paper for reference).  One runs
% this program by putting the folder containing this file AND the file
% 'fingIPVeqns.m' in the Matlab's directory and typing 'fingIPVint(1)' into
% the command window.

global p m a b c f g h 

%par=[1/45 144 36 12 12];

%known constants and variable paramters
m=1/45;
a=.3;
b=144;
c=12;
f=0.5;
g=1;
h=.2;
p=0.6875;
col='';


x0=[.5,.25,.25,0,0]; %initial conditions

s0=[];
s1=[];
w0=[];
w1=[];
r=[];
T=[];




op=odeset('reltol',1e-10, 'abstol',1e-11,'stats','on');

tic

   tspan=[0, tf];
    [t,x] = ode45('fingIPVeqns',tspan,x0,op);
     s0=[s0;x(:,1)];
     s1=[s1;x(:,2)];
     w0=[w0;x(:,3)];
     w1=[w1;x(:,4)];
     r=[r;x(:,5)];
     T=[T;t];
toc 
   % figure
    plot(T,s0,T,s1,T,w0,T,w1)
    legend('sus','sus vac','wild','wild vac','immune',0);
    title(['optimal vaccinations => (' num2str(x(end,1)) ', ' num2str(x(end,2)) ',' num2str(x(end,3)) ', ' num2str(x(end,4)) ')']);
    xlabel('time');
    ylabel('% population');    
   % hold off;
    final=[s0(end) s1(end) w0(end) w1(end) r(end)]