function harvest
clc,clf, close  all

global q r K E p c x0 T
q=10^(-5);
r=0.05;
K=250000;
E=5000;
p=1;
c=1;
x0=150000;
T=10; % final time
PsiTrange =[];% keeps the terminal values for psi(T)
Psi0range = -0.99:0.01:0.99;

hold on


for psi0value = Psi0range

    [t,X,U,J,psiTvalue] = shoot(psi0value);
    
    PsiTrange=[PsiTrange;psiTvalue];
   
    plot(X(:,1),X(:,2),'-')
        axis([0 K -1 1])
        grid on
        xlabel('x'), ylabel('psi')
        title('Phase portrait with switching function')
        %text(1.2, 0.1, ['tfinal = ' num2str(tfinal)])
end
    hold off
    
    figure(2)
    plot(Psi0range,PsiTrange,'-o')
    grid on, xlabel('Initial Psi(0)'), ylabel('Final Psi(T)')
    title('Result of the shooting method')
    
    [minPsi,ind] = min(abs(PsiTrange)); % finds the value of PsiT closest to 0
    
    psi0 = Psi0range(ind);
    fprintf('Approximate value of optimal psi0 is %g\n\n',psi0)
    
    [t,X,U,J,psiT] = shoot(psi0);
    
    figure(3)
    subplot(2,1,1)
    plot(t,U,'-*');
    xlabel('t'), ylabel('u*')
    title('Optimal Fishing effort u=u*(t)')
    subplot(2,1,2)
    plot(t,X(:,1),'-*');
    xlabel('t'), ylabel('x*')
    title('Resulting Fish poluation levels x=x*(t)')
    fprintf('The Optimal Profit is %6.2f dollars\n', J)
    fprintf('See figure for the optimal fishing effort and resulting fish polulation') 
end
    
function [t,X,U,J,PsiT] = shoot(Psi0)

global q r K E p c x0 T

            X0 = [x0; Psi0];
            Ts=[]; % keeps the switching times (if any)
            if Switch(0,X0)<0 
                u=0;
            else u=E;
            end

            RHS=@(t,X) [r*X(1)*(1-X(1)/K)-q*u*X(1);
                        -(p*q*u+(r*(1-2*X(1)/K)-q*u)*X(2))];

            opts = odeset('events',@Switch);

            [t,X,ts] = ode45(RHS,[0 T],X0,opts);
            % ts is the switch time
            U=u*ones(size(t));

            if size(ts)>0

                u=E-u; 
                % changes the value of u from 0 to E or viceversa.
                RHS=@(t,X) [r*X(1)*(1-X(1)/K)-q*u*X(1);
                        -(p*q*u+(r*(1-2*X(1)/K)-q*u)*X(2))];

                opts = odeset('events',@Switch);

                [tnew,Xnew,tsnew] = ode45(RHS,[ts T],X(end,:),opts);

                Ts=[Ts;tsnew];
                ts=tsnew;

                t=[t;tnew];
                X=[X;Xnew];
                U=[U; u*ones(size(tnew))];

            end
            J=trapz(p*q*U.*X(:,1)-c*U);
            PsiT=X(end,2);          
end


function [Sstop,isterminal,direction] = Switch(t,X)
 
global q p c 
    Sstop = q*X(1)*(p-X(2))-c;
    isterminal = 1;
    direction = 0; % records all zeros
end