function [T,Y]=mainsci

%this is the main function for solving the OPV model optimally.  
% To use this program, simply put this program AS WELL AS shootsci.m in the
% file directory, and type 'mainsci' into Matlab's command window.

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

beta=[0 0 0 1];  %desired final values for costate variables
tolerance=.18; %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 1 0 0 0 0 1 0 0 0 0 1];
o0=[1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1];
%put the initial guess here  (quess where you want to end)
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('shootsci',[0 2],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));
  %      axis([0 .8 0 0.02])
        legend('sus','vacc','wild','rem',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,5:8)';   %the final values of the costate variables
    Phi=q-beta';
    Phiprime=[Y(r,25:28); Y(r,29:32); Y(r,33:36); Y(r,37:40)];
    znew=z-Phiprime\Phi;
    error=norm(znew-z)
    z=znew;
 toc
    [a,b]=size(Y(:,5));   
    sw=par(1)*(Y(:,5)-Y(:,6));  %let's plot the control values!!
       for i=1:1:a
           if sw(i) < 0 
               p=.6875;
           else if sw(i) > 0 
               p=0;
           else p=.5;
               end;
           end;
            P=[P; p];
        end
   plot(T,P)
end;