% First_Order_System_2

% Analysis of system of autonomous first-order ODE's:
%   dx/dt = F(x,y)
%   dy/dt = G(x,y)

% Uses Euler algorithm and/or one of the Matlab solvers

%  APMA 0410
%  Fall 2012

clear

% CHOOSE FROM THE FOLLOWING (AT MOST ONE MATLAB SOLVER):
Euler_flag  = 1;     % set to 1 to use Euler scheme
ode45_flag  = 0;     % set to 1 to use ode45
ode23_flag  = 0;     % set to 1 to use ode23
ode15s_flag = 0;     % set to 1 to use ode15s

if ode45_flag + ode23_flag + ode15s_flag > 1
    error('Use at most one Matlab solver')
end
                        
dt = .001;     % time step for Euler's method

% DEFINE the SYSTEM
%    CAUTION: FOR SOME OF THE SYSTEMS, BOUNDING BOX MAY HAVE TO BE ADJUSTED

% linear system - center - note the inaccuracy of Euler, but also of Matlab solvers
% F = @(x,y)    x + 3*y ; 
% G = @(x,y) -2*x -   y ; 
% min_x = -5;
% max_x =  5;
% min_y = -5;
% max_y =  5;
% t_max = 10;     % length of simulation

% Example 1 from Lecture Notes (Non-Linear Systems of Autonomous First-Order DE's -- Examples)
% F = @(x,y) y ; 
% G = @(x,y) x.*(x - 1) ;
% min_x = -0.8;
% max_x =  1.8;
% min_y = -0.8;
% max_y =  0.8;
% min_x = -500;
% max_x =  500;
% min_y = -500;
% max_y =  500;
% t_max = 20;     % length of simulation


% Example 2 from Lecture Notes (Non-Linear Systems of Autonomous First-Order DE's -- Examples)
% van der Pol oscillator
% mu = 1;
% F = @(x,y) y ; 
% G = @(x,y) mu*(1 - x.^2).*y - x ; 
% min_x = -5;
% max_x =  5;
% min_y = min(-5,-2*mu);
% max_y = max( 5, 2*mu);
% t_max = 10*mu;     % length of simulation
% Remark: for large mu (e.g. 100), the van der Pol oscillator is stiff:
%     Euler's method badly overshoots the slow "relaxation branches"
%        of the attractor when started far from the origin,
%     ODE45 takes forever,
%     ODE15s works superbly (set t_max=1000).

% Example 3 from Lecture Notes (Non-Linear Systems of Autonomous First-Order DE's -- Examples)
% Do not use ODE45 or ODE23 for large t_max (e.g. 20)
% F = @(x,y) x.*y ; 
% G = @(x,y) y - x ; 
% min_x = -1;
% max_x =  2;
% min_y = -1;
% max_y =  2;
% t_max = 10;     % length of simulation


% Example 4 from Lecture Notes (Non-Linear Systems of Autonomous First-Order DE's -- Examples)
% Parabolic V-nullcline - No equilibrium
% F = @(x,y) x.^2 - y + 2; 
% G = @(x,y) y + 2*x ;
% min_x = -4;
% max_x =  2;
% min_y = -1;
% max_y =  8;
% t_max = 10;     % length of simulation


% Example 5 from Lecture Notes (Non-Linear Systems of Autonomous First-Order DE's -- Examples)
% F = @(x,y) 1 + y + x.^2 ; 
% G = @(x,y) y/2 - 3*x ; 
% min_x = -10;
% max_x =  10;
% min_y = -45;
% max_y =   5;
% t_max = 20;     % length of simulation

% Reverse the time:
% F = @(x,y) -(1 + y + x.^2) ; 
% G = @(x,y) -(y/2 - 3*x) ; 
% min_x = -5;
% max_x =  5;
% min_y = -10;
% max_y =  5;
% t_max = 10;     % length of simulation


% Saddle-node bifurcation at mu = 0 
% mu = 1;
% F = @(x,y) x.*x - mu ; 
% G = @(x,y) - y ; 
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 10;     % length of simulation


% saddle and two spiral sinks
% F = @(x,y) y ; 
% G = @(x,y) - y + x - x.^3 ; 
% min_x = -5;
% max_x =  5;
% min_y = -10;
% max_y =  10;
% t_max = 20;     % length of simulation 

% species-competition model
% F = @(x,y) (1.0 - 1.0*x - 1.0*y).*x; 
% G = @(x,y) (4.0 - 7.0*x - 3.0*y).*y;  
% min_x = -0.2;
% max_x =  1.5;
% min_y = -0.2;
% max_y =  1.5;
% t_max = 10;     % length of simulation 


% Lotka-Volterra (don't click on negative values)
% F = @(x,y) x.*( 2.0 - 1.0*y) ; 
% G = @(x,y) y.*(-3.0 + 0.5*x) ; 
% min_x = -0.2;
% max_x =  10;
% min_y = -0.2;
% max_y =  10;
% t_max = 10;     % length of simulation 


% Lotka-Volterra with harvesting:
% F = @(x,y) x.*( 2.0 - 1.0*y) - 0.5*x ; 
% G = @(x,y) y.*(-3.0 + 0.5*x) - 0.6*y ; 
% min_x = -0.2;
% max_x =  10;
% min_y = -0.2;
% max_y =  3;
% t_max = 10;     % length of simulation 


% Lotka-Volterra with logistic terms:
% F = @(x,y) x.*( 3.0 - 2.0*y) - 0.2*x.^2; 
% G = @(x,y) y.*(-3.0 + 2.0*x) - 0.2*y.^2; 
% min_x = -0.2;
% max_x =  10;
% min_y = -0.2;
% max_y =  3;
% t_max = 10;     % length of simulation 


% FitzHugh-Nagumo Model
% Input = .7; %        Hopf bifurcations around .34 and around 1.42
% F = @(x,y) x - x.^3/3 - y + Input ;
% G = @(x,y) 0.08.*(x + 0.7 - 0.8.*y) ;
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 100;     % length of simulation 


% F = @(x,y) x - x.^3/3 - y ;
% G = @(x,y) x  ;
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 10;     % length of simulation 

% Limit cycle:
% F = @(x,y) y ;
% G = @(x,y) y - y.^3/3 - x ;
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 20;     % length of simulation 


% F = @(x,y)  y - .1*(x.^2 + y.^2) .*x ; 
% G = @(x,y) -x - .1*(x.^2 + y.^2) .*y ; 
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 100;     % length of simulation 


% Harmonic oscillator
% omega = 2;
% F = @(x,y) - omega.*y; 
% G = @(x,y)   omega.*x; 
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 10;     % length of simulation 
% equivalent to:
%   r'     = 0
%   theta' = omega 
% where:
%  x = r*cos(theta)
%  y = r*sin(theta)


% Illustration of Poincare-Bendixson Theorem
% omega = 2;
% F = @(x,y) x - omega.*y - x.*(x.^2 + y.^2); 
% G = @(x,y) omega.*x + y - y.*(x.^2 + y.^2); 
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 10;     % length of simulation 
% equivalent to:
%   r'     = r - r^3
%   theta' = omega 
% where:
%  x = r*cos(theta)
%  y = r*sin(theta)

% Supercritical Hopf bifurcation at mu = 0
% mu = .1;
% omega = 4;
% F = @(x,y)  mu.*x - omega.*y - x.*(x.^2 + y.^2); 
% G = @(x,y)  omega.*x + mu.*y - y.*(x.^2 + y.^2); 
% min_x = -2;
% max_x =  2;
% min_y = -2;
% max_y =  2;
% t_max = 100;     % length of simulation 
% equivalent to:
%   r'     = mu*r - r^3
%   theta' = omega 
% where:
%  x = r*cos(theta)
%  y = r*sin(theta)

% Subcritical Hopf bifurcation at mu = 0
% mu = -1;
% omega = 4;
% F = @(x,y)  mu.*x - omega.*y + x.*(x.^2 + y.^2); 
% G = @(x,y)  omega.*x + mu.*y + y.*(x.^2 + y.^2); 
% min_x = -2;
% max_x =  2;
% min_y = -2;
% max_y =  2;
% t_max = 100;     % length of simulation 
% equivalent to:
%   r'     = mu*r + r^3
%   theta' = omega 
% where:
%  x = r*cos(theta)
%  y = r*sin(theta)


% Excitatory-Inhibitory System (from EB+DL 1999) 
% F = @(x,y) - x + .5*tanh(10*x - 10*y); 
% G = @(x,y) - y + .5*tanh(10*x - 1*y); 
% min_x = -1;
% max_x =  1;
% min_y = -1;
% max_y =  1;
% t_max = 10;     % length of simulation 
 
% 
% F = @(x,y) y.^2 - x - 2;
% G = @(x,y) y - y.^3/3 - x ;
% min_x = -3;
% max_x =  3;
% min_y = -3;
% max_y =  3;
% t_max = 20;     % length of simulation 


% HW 6.4
% F = @(x,y) (x + 1).*y;
% G = @(x,y) y - y.^2 - x + 1;
% min_x = -30;
% max_x =  30;
% min_y = -30;
% max_y =  30;
% t_max = 20;     % length of simulation 

% nice for final
F = @(x,y)   y.^2 - x - 2;
G = @(x,y) - y.^2 - x ;
min_x = -5;
max_x =  5;
min_y = -5;
max_y =  5;
t_max = 20;     % length of simulation 

% nice for final
% F = @(x,y) y.^2 - x.^2;
% G = @(x,y)  -y.^2 - x + 2;
% min_x = -5;
% max_x =  5;
% min_y = -5;
% max_y =  5;
% t_max = 20;     % length of simulation 


% nice for final
% F = @(x,y) y.^2 - x.^2;
% G = @(x,y)  -y.^2 - x.^2 + 2;
% min_x = -20;
% max_x =  20;
% min_y = -20;
% max_y =  20;
% t_max = 20;     % length of simulation 



F_and_G = @(t, Y) [F(Y(1),Y(2)) ; G(Y(1),Y(2))];


% DIRECTION FIELD 

% Define the grid points:
[x,y] = meshgrid(linspace(min_x,max_x,15),linspace(min_y,max_y,15));

figure(1), set(gca,'fontsize',16)
% The derivatives at all grid points:
dx = F(x,y);
dy = G(x,y);

% Put the (dx,dy) arrow at each grid point (x,y):
quiver(x,y,dx,dy)
axis([min_x max_x min_y max_y])
grid on
axis image
set(gca,'fontsize',16)
xlabel('x')
ylabel('y')
hold on

% plot the nullclines:
h_V_nc = ezplot(F,[min_x,max_x,min_y,max_y]);   % vertical   NC
h_H_nc = ezplot(G,[min_x,max_x,min_y,max_y]);   % horizontal NC
set(h_V_nc,'color','g','linewidth',2)
set(h_H_nc,'color','m','linewidth',2)

% TRAJECTORIES

% To create a new trajectory, click mouse on desired initial point in (x,y) plane.

% To end simulation (e.g. to change parameter values),
%   click mouse in grey margin, i.e., outside of permissible area.


time_axis = 0:dt:t_max;
n_time_steps = length(time_axis);

X = zeros(2,n_time_steps);
X(1,:) = .5*(min_x + max_x)*ones(1,n_time_steps);
X(2,:) = .5*(min_y + max_y)*ones(1,n_time_steps);
% this array will contain an entire trajectory

h_E = [];

while X(1,1)>=min_x && X(1,1)<=max_x && X(2,1)>=min_y && X(2,1)<=max_y  
    % will create a new trajectory until mouse is clicked in margin

    figure(1), set(gca,'fontsize',16)
    X(:,1) = ginput(1);    % click mouse on desired initial weight vector

    if Euler_flag       % Do Euler's method
        
        i = 1;
        while i <= n_time_steps && X(1,i)>10*min_x && X(1,i)<10*max_x && X(2,i)>10*min_y && X(2,i)<10*max_y
            i = i + 1;
            X(:,i) = X(:,i-1) + [F(X(1,i-1),X(2,i-1));G(X(1,i-1),X(2,i-1))]*dt;
        end

        tt = 1:i-1;
        t_a = time_axis(tt);

        if ~ isempty(tt)

            h_E = plot(X(1,tt),X(2,tt),'r','linewidth',2);
            figure(2)
            plot(t_a,X(1,tt),t_a,X(2,tt),'linewidth',2)
            legend('x(t)','y(t)')
            xlabel('t')
            grid

        end

    end

    if ode45_flag
        [t,Y] = ode45(F_and_G,[0 t_max],X(:,1));
    end
    if ode23_flag
        [t,Y] = ode23(F_and_G,[0 t_max],X(:,1));
    end
    if ode15s_flag
        [t,Y] = ode15s(F_and_G,[0 t_max],X(:,1));
    end

    h_M = [];
    if ode45_flag || ode23_flag || ode15s_flag
        X_m = Y(:,1);
        Y_m = Y(:,2);
        figure(1), set(gca,'fontsize',16)
        h_M = plot(X_m,Y_m,'k','linewidth',1);
    end
    
end

figure(1), set(gca,'fontsize',16)
x_deriv = func2str(F);
x_deriv_name = x_deriv(7:end);
y_deriv = func2str(G);
y_deriv_name = y_deriv(7:end);
title(['x'' = ' x_deriv_name '       y'' = ' y_deriv_name])
if isempty(h_E) && ~ isempty(h_M) && ~ isempty(h_V_nc) && ~ isempty(h_H_nc)
legend([h_M h_V_nc(1) h_H_nc(1)],'Matlab','V-nc','H-nc','Location','Best')
end
if isempty(h_M) && ~ isempty(h_E) && ~ isempty(h_V_nc) && ~ isempty(h_H_nc)
legend([h_E h_V_nc(1) h_H_nc(1)],'Euler','V-nc','H-nc','Location','Best')
end
if ~isempty(h_M) && ~ isempty(h_E) && ~ isempty(h_V_nc) && ~ isempty(h_H_nc)
legend([h_E h_M h_V_nc(1) h_H_nc(1)],'Euler','Matlab','V-nc','H-nc','Location','Best')
end
hold off

