function [T,Y]=mainIPV

% This is the main function for solving the IPV model optimally.  One must
% add the folder containing this AS WELL AS the file 'shootIPV.m' in the
% Matlab's file directory.  One runs this program by typing 'mainIPV' into
% Matlab's command window.

par=[1/45 .3 144 12 0.5 1 .2 10 10];
x0=[.75 0 .25 0 0];

beta=[0 0 1 1 1];  %desired final values for costate variables
tolerance=5e-005; %error tolerance for Newton's method
error=1;    %make the error large so that we do at least one iteration of Newton's method

n0=[1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1];
o0=[1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1];
%put the initial guess here  (quess where you want to end up!)
z=beta';  



while (error > tolerance)  %newton's loop
 tic  %time it all
    P=[];  %definitions for later
    i=0;
    xinit=[x0 z'  n0 o0];
    [T,Y]=ode45('shootIPV',[0 1],xinit,[],par); %the non-stiff package is more efficient
        figure, hold on
        plot(T,Y(:,1),T,Y(:,2),T,Y(:,3),T,Y(:,4),T,Y(:,5));
%        axis([0 .8 0 0.02])
        legend('sus','sus vac','wild vac','wild','immune',0);
        title(['optimal vaccinations => (' num2str(Y(end,1)) ', ' num2str(Y(end,2)) ',' num2str(Y(end,3)) ', ' num2str(Y(end,4)) ')']);
        xlabel('time');
        ylabel('% population');
        drawnow
    [r,c]=size(Y);
    q=Y(r,6:10)';   %the final values of the costate variables
    Phi=q-beta';
    Phiprime=[Y(r,36:40); Y(r,41:45); Y(r,46:50); Y(r,51:55); Y(r,56:60)];
    znew=z-Phiprime\Phi;
    error=norm(znew-z)
    z=znew;
 toc
    [a,b]=size(Y(:,6));   
    sw=par(1)*Y(:,6)+par(1)*par(2)*(Y(:,7)-Y(:,10));  %let's plot the control values!!
       for i=1:1:a
           if sw(i) < 0 
               p=.9;
           else if sw(i) > 0 
               p=0;
           else p=.5;
               end;
           end;
            P=[P; p];
        end
   plot(T,P)
end;