close

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

A = input('Enter a value for the weight parameter A: ');
disp('         ')
B = input('Enter a value for the weight parameter B: ');
disp('        ')
while(flag1==0)
    var1 = input('Enter a strictly negative value for the weight parameter C: ');
    if(var1<0)
        C = var1;
        flag1=1;
        disp('        ')
    else
        disp('        ')
        disp('ERROR: C must be negative')
        disp('     ')
    end
end
flag=0;
x0 = input('Enter the initial condition of the state (x_0): ');
disp('           ')
while(flag2==0)
    var2 = input('Enter the upper bound M: ');
    if(var2>0)
        M = var2;
        flag2=1;
        disp('           ')
    else
        disp('       ')
        disp('ERROR: M must be positive')
        disp('      ')
    end
end

h=2/500;
t = linspace(0,2,501);
a1 = A*(exp(2-t)-1);
p = (-a1-B)/(2*C);
u1 = min(M,max(0,p));

s1 = zeros(1,501);
s1(1) = x0;
for i=1:500
    s1(i+1) = s1(i) + (s1(i)+u1(i))*h;
end

while(flag3==0)
    disp('Would you like to plot')
    disp('   1. state on optimal path (x^*)')
    disp('   2. adjoint (lambda)')
    disp('   3. optimal control (u^*)')
    disp('   4. state and control')
    disp('   5. state and adjoint')
    disp('   6. control and adjoint')
    disp('   7. all three')
    var3=input('Type 1, 2, 3, 4, 5, 6, or 7: ');
    
    if((var3==1) | (var3==2) | (var3==3) | (var3==4) | (var3==5) | (var3==6) | (var3==7)) 
        disp('       ')
        flag3=1;
        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)
                flag4=1;
                disp('          ')
                while(flag5==0)
                    disp('Which parameter would you like to vary?')
                    disp('   1. A')
                    disp('   2. B')
                    disp('   3. C')
                    disp('   4. x_0')
                    disp('   5. M')
                    var5=input('Type 1, 2, 3, 4, or 5: ');
                    if(var5==1)
                        flag5=1;
                        disp('     ')
                        A2=input('Enter a second A value: ');
                        a2 = A2*(exp(2-t)-1);
                        p2 = (-a2-B)/(2*C);
                        u2 = min(M,max(0,p2));
                        s2 = zeros(1,501);
                        s2(1) = x0;
                        for i=1:500
                            s2(i+1) = s2(i) + (s2(i)+u2(i))*h;
                        end
                    elseif(var5==2)
                        flag5=1;
                        disp('       ')
                        B2=input('Enter a second B value: ');
                        a2 = A*(exp(2-t)-1);
                        p2 = (-a2-B2)/(2*C);
                        u2 = min(M,max(0,p2));
                        s2 = zeros(1,501);
                        s2(1) = x0;
                        for i=1:500
                            s2(i+1) = s2(i) + (s2(i)+u2(i))*h;
                        end
                    elseif(var5==3)
                        flag5=1;
                        disp('        ')
                        while(tempflag==0)
                            temp=input('Enter a second C value: ');
                            if(temp<0)
                                C2=temp;
                                tempflag=1;
                                a2 = A*(exp(2-t)-1);
                                p2 = (-a2-B)/(2*C2);
                                u2 = min(M,max(0,p2));
                                s2 = zeros(1,501);
                                s2(1) = x0;
                                for i=1:500
                                    s2(i+1) = s2(i) + (s2(i)+u2(i))*h;
                                end
                            else
                                disp('         ')
                                disp('ERROR: C must be negative')
                                disp('         ')
                            end
                        end
                    elseif(var5==4)
                        flag5=1;
                        disp('        ')
                        x02=input('Enter a second x_0 value: ');
                        a2 = A*(exp(2-t)-1);
                        p2 = (-a2-B)/(2*C);
                        u2 = min(M,max(0,p2));
                        s2 = zeros(1,501);
                        s2(1) = x02;
                        for i=1:500
                            s2(i+1) = s2(i) + (s2(i)+u2(i))*h;
                        end
                    elseif(var5==5)
                        flag5=1;
                        disp('          ')
                        while(tempflag2==0)
                            temp2=input('Enter a second M value: ');
                            if(temp2>0)
                                M2=temp2
                                tempflag2=1;
                                a2 = A*(exp(2-t)-1);
                                p2 = (-a2-B)/(2*C);
                                u2 = min(M2,max(0,p2));
                                s2 = zeros(1,501);
                                s2(1) = x0;
                                for i=1:500
                                    s2(i+1) = s2(i) + (s2(i)+u2(i))*h;
                                end
                            else
                                disp('          ')
                                disp('ERROR: M must be positive')
                                disp('         ')
                            end
                        end
                        
                    else
                        disp('        ')
                        disp('Pardon?')
                        disp('        ')
                        
                    end
                end 
                
            elseif(var4==2)
                flag4=1;
                if(var3==1)
                    plot(t,s1)
                    xlabel('Time')
                    ylabel('State')
                elseif(var3==2)
                    plot(t,a1)
                    xlabel('Time')
                    ylabel('Adjoint')
                elseif(var3==3)
                    plot(t,u1)
                    xlabel('Time')
                    ylabel('Control')
                    axis([0 2 -1 M+1])
                elseif(var3==4)
                    subplot(2,1,1);plot(t,s1)
                    subplot(2,1,1);xlabel('Time')
                    subplot(2,1,1);ylabel('State')
                    subplot(2,1,2);plot(t,u1)
                    subplot(2,1,2);xlabel('Time')
                    subplot(2,1,2);ylabel('Control')
                    subplot(2,1,2);axis([0 2 -1 M+1])
                elseif(var3==5)
                    subplot(2,1,1);plot(t,s1)
                    subplot(2,1,1);xlabel('Time')
                    subplot(2,1,1);ylabel('State')
                    subplot(2,1,2);plot(t,a1)
                    subplot(2,1,2);xlabel('Time')
                    subplot(2,1,2);ylabel('Adjoint')
                elseif(var3==6)
                    subplot(2,1,1);plot(t,u1)
                    subplot(2,1,1);xlabel('Time')
                    subplot(2,1,1);ylabel('Control')
                    subplot(2,1,1);axis([0 2 -1 M+1])
                    subplot(2,1,2);plot(t,a1)
                    subplot(2,1,2);xlabel('Time')
                    subplot(2,1,2);ylabel('Adjoint')
                elseif(var3==7)
                    subplot(3,1,1);plot(t,s1)
                    subplot(3,1,1);xlabel('Time')
                    subplot(3,1,1);ylabel('State')
                    subplot(3,1,2);plot(t,u1)
                    subplot(3,1,2);xlabel('Time')
                    subplot(3,1,2);ylabel('Control')
                    subplot(3,1,2);axis([0 2 -1 M+1])
                    subplot(3,1,3);plot(t,a1)
                    subplot(3,1,3);xlabel('Time')
                    subplot(3,1,3);ylabel('Adjoint')
                end    
             
            else
                disp('         ')
                disp('Pardon?')
                disp('          ')
                
            end    
            
        end
        
    else
        disp('         ')
        disp('Pardon?')
        disp('         ')
        
    end
    
end        
            
if(var4==1)
    if(var3==1)
        plot(t,s1,'b',t,s2,'g')
        xlabel('Time')
        ylabel('States')
        legend('First Value','Second Value',0)
    elseif(var3==2)
        plot(t,a1,'b',t,a2,'g')
        xlabel('Time')
        ylabel('Adjoints')
        legend('First Value','Second Value',0)
    elseif(var3==3)
        plot(t,u1,'b',t,u2,'g')
        xlabel('Time')
        ylabel('Controls')
        legend('First Value','Second Value',0)
        if(var5==5)
            axis([0 2 -1 max(M+1,M2+1)])
        else
            axis([0 2 -1 M+1])
        end
    elseif(var3==4)
        subplot(2,1,1);plot(t,s1,'b',t,s2,'g')
        subplot(2,1,1);xlabel('Time')
        subplot(2,1,1);ylabel('States')
        subplot(2,1,1);legend('First Value','Second Value',0)
        subplot(2,1,2);plot(t,u1,'b',t,u2,'g')
        subplot(2,1,2);xlabel('Time')
        subplot(2,1,2);ylabel('Controls')
        subplot(2,1,2);legend('First Value','Second Value',0)
        if(var5==5)
            subplot(2,1,2);axis([0 2 -1 max(M+1,M2+1)])
        else
            subplot(2,1,2);axis([0 2 -1 M+1])
        end
    elseif(var3==5)
        subplot(2,1,1);plot(t,s1,'b',t,s2,'g')
        subplot(2,1,1);xlabel('Time')
        subplot(2,1,1);ylabel('States')
        subplot(2,1,1);legend('First Value','Second Value',0)
        subplot(2,1,2);plot(t,a1,'b',t,a2,'g')
        subplot(2,1,2);xlabel('Time')
        subplot(2,1,2);ylabel('Adjoints')
        subplot(2,1,2);legend('First Value','Second Value',0)
    elseif(var3==6)
        subplot(2,1,1);plot(t,u1,'b',t,u2,'g')
        subplot(2,1,1);xlabel('Time')
        subplot(2,1,1);ylabel('Controls')
        subplot(2,1,1);legend('First Value','Second Value',0)
        if(var5==5)
            subplot(2,1,1);axis([0 2 -1 max(M+1,M2+1)])
        else
            subplot(2,1,1);axis([0 2 -1 M+1])
        end
        subplot(2,1,2);plot(t,a1,'b',t,a2,'g')
        subplot(2,1,2);xlabel('Time')
        subplot(2,1,2);ylabel('Adjoints')
        subplot(2,1,2);legend('First Value','Second Value',0)
    elseif(var3==7)
        subplot(3,1,1);plot(t,s1,'b',t,s2,'g')
        subplot(3,1,1);xlabel('Time')
        subplot(3,1,1);ylabel('States')
        subplot(3,1,1);legend('First Value','Second Value',0)
        subplot(3,1,2);plot(t,u1,'b',t,u2,'g')
        subplot(3,1,2);xlabel('Time')
        subplot(3,1,2);ylabel('Controls')
        subplot(3,1,2);legend('First Value','Second Value',0)
        if(var5==5)
            subplot(3,1,2);axis([0 2 -1 max(M+1,M2+1)])
        else
            subplot(3,1,2);axis([0 2 -1 M+1])
        end
        subplot(3,1,3);plot(t,a1,'b',t,a2,'g')
        subplot(3,1,3);xlabel('Time')
        subplot(3,1,3);ylabel('Adjoints')
        subplot(3,1,3);legend('First Value','Second Value',0)
    end
end
