close

flag1=0;
flag2=0;
flag3=0;
flag4=0;
flag5=0;

while(flag1==0)
    var1 = input('Enter the initial number of fish (x_0): ');
    if(var1<=0)
        disp('           ')
        disp('ERROR: x_0 must be positive.')
        disp('            ')
    else
        flag1=1;
        x0=var1;
    end    
end
disp('           ')
while(flag2==0)
    var2 = input('Enter a value for the time (T): ');
    if(var2<=0)
        disp('           ')
        disp('ERROR: T must be positive.')
        disp('            ')
    else
        flag2=1;
        T=var2;
    end    
end

disp('               ')

t=linspace(0,T,500);
a=2*tan(-2*t + 2*T);
x1=(cos(2*T - 2*t)).^2;
k=x1(1);
c=x0/k;
x=c*x1;
u=x.*(2-a).^2;
r=2*(u.^.5).*(x.^.5);

while(flag3==0)
    disp('Would you like to plot')
    disp('   1. Number of fish (state)')
    disp('   2. Intensity of harvesting (control)')
    disp('   3. Harvesting rate (r)')
    disp('   4. all three')
    var3 = input('Type 1, 2, 3, or 4: ');

    if((var3==1) | (var3==2) | (var3==3) | (var3==4))
        flag3=1;
        disp('         ')
        while(flag4==0)
            disp('Would you like to vary any parameters?')
            disp('1. Yes')
            disp('2. No')
            var4 = input('Type 1 or 2: ');
            
            if(var4==1)
                disp('         ')
                flag4 = 1;
                while(flag5==0)
                    disp('Which parameter would you like to vary?')
                    disp('1. x_0')
                    disp('2. T')
                    var5=input('Type 1 or 2: ');
                    if(var5==1)
                        disp('        ')
                        x02=input('Enter a second x_0 value: ');
                        t2=t;
                        c2=x02/k;
                        x2=c2*x1;
                        u2=x2.*(2-a).^2;
                        r2=2*(u2.^.5).*(x2.^.5);
                        flag5=1;
                    elseif(var5==2)
                        disp('        ')
                        T2=input('Enter a second T value: ');
                        t2=linspace(0,T2,500);
                        a2=2*tan(-2*t2 + 2*T2);
                        x12=(cos(2*T2-2*t2)).^2;
                        k2=x12(1);
                        c2=x0/k2;
                        x2=c2*x12;
                        u2=x2.*(2-a2).^2;
                        r2=2*(u2.^.5).*(x2.^.5);
                        flag5=1;
                    else
                        disp('         ')
                        disp('Pardon?')
                        disp('           ')
                    end
                end
            
            elseif(var4==2)
                disp('            ')
                flag4=1;
                if(var3==1)
                    plot(t,x)
                    xlabel('Time')
                    ylabel('Fish Population')
                elseif(var3==2)
                    plot(t,u)
                    xlabel('Time')
                    ylabel('Harvesting Intensity')
                elseif(var3==3)
                    plot(t,r)
                    xlabel('Time')
                    ylabel('Harvesting Rate')
                elseif(var3==4)
                    subplot(3,1,1);plot(t,x)
                    subplot(3,1,1);xlabel('Time')
                    subplot(3,1,1);ylabel('Fish')
                    subplot(3,1,2);plot(t,u) 
                    subplot(3,1,2);xlabel('Time')
                    subplot(3,1,2);ylabel('Harvesting Intensity')
                    subplot(3,1,3);plot(t,r)
                    subplot(3,1,3);xlabel('Time')
                    subplot(3,1,3);ylabel('Harvesting Rate')
                end
                
            else
                disp('           ')
                disp('Pardon?')
                disp('            ')
            end
        end
        
    else 
        disp('     ')
        disp('Pardon?')
        disp('          ')
    end
end

if(var4==1)
    if(var3==1)
        plot(t,x,'b',t2,x2,'g')
        xlabel('Time')
        ylabel('Fish Population')
        legend('First value','Second value',0)
    elseif(var3==2)
        plot(t,u,'b',t2,u2,'g')
        xlabel('Time')
        ylabel('Harvesting Intensity')
        legend('First value','Second value',0)
    elseif(var3==3)
        plot(t,r,'b',t2,r2,'g')
        xlabel('Time')
        ylabel('Harvesting Rate')
        legend('First value','Second value',0)
    elseif(var3==4)
        subplot(3,1,1);plot(t,x,'b',t2,x2,'g')
        subplot(3,1,1);legend('First value','Second value',0)
        subplot(3,1,1);ylabel('Fish')
        subplot(3,1,1);xlabel('Time')
        subplot(3,1,2);plot(t,u,'b',t2,u2,'g')
        subplot(3,1,2);xlabel('Time')
        subplot(3,1,2);ylabel('Harvesting Intensity')
        subplot(3,1,2);legend('First value','Second value',0)
        subplot(3,1,3);plot(t,r,'b',t2,r2,'g')
        subplot(3,1,3);xlabel('Time')
        subplot(3,1,3);ylabel('Harvesting Rate')
        subplot(3,1,3);legend('First value','Second value',0)
    end
end
