% Linear_First_Order_System.m

% Analysis of system of autonomous linear first-order ODE's:

%   dx/dt = a*x + b*y + e
%   dy/dt = c*x + d*y + f

%  APMA 0410
%  Fall 2012

clear

% PARAMETERS

% a =  1; b =  0; c =  0; d =  1;
% a =  1; b =  0; c =  0; d =  2;
% a = -1; b =  0; c =  0; d = -2;
% a =  1; b = -2; c =  1; d =  1;
% a = -1; b =  2; c = -1; d = -1;
% a =  1; b =  0; c =  0; d = -1;
% a =  1; b =  0; c =  0; d = -5;
% a =  2; b =  2; c =  1; d =  3;
% a =  1; b =  3; c = -2; d = -1;
% a =  0; b =  0; c =  0; d =  1;
% a =  1; b =  1; c =  1; d =  1;
% a =  3; b =  1; c =  0; d =  3;
% a = 11; b = -9; c =  1; d =  5;
% a =  1; b = -2; c =  2; d = -4;
% a =  3; b = -2; c =  2; d = -1;
% a =  1; b = -2; c = .5; d = -1;
% a =  3/4; b =  1; c =  3/4-49/64; d =  1;
% a =  2; b =  2; c =  3; d =  1;
% a = -2; b = -2; c = -1; d = -3;
% a = -1; b =  0; c = -1; d =  1;
% a =  1; b = -3; c =  2; d = -3;
a =  1; b =  1; c =  0; d =  2;

e =  0; f =  0;
% e =  4; f =  0;
 
F = @(x,y) a*x + b*y + e ; 
G = @(x,y) c*x + d*y + f ; 

% Matrix form:
A = [a b
     c d];
non_homogeneous_term = [e;f];
 

%% Eigenvalue Analysis

[Eigenvectors,Eigenvalues] = eig(A);

eigenvalue_1 = Eigenvalues(1,1);
eigenvalue_2 = Eigenvalues(2,2);

eigenvector_1 = Eigenvectors(:,1);
eigenvector_2 = Eigenvectors(:,2);

disp(['First Eigenvalue : ',num2str(eigenvalue_1)])
disp(['First Eigenvector : ',num2str(eigenvector_1')])
disp(' ')
disp(['Second Eigenvalue : ',num2str(eigenvalue_2)])
disp(['Second Eigenvector : ',num2str(eigenvector_2')])
    

%% Trace-Determinant Analysis

T = a + d;       % trace of matrix A
D = a*d - b*c;   % determinant of matrix A

% Visualize the parameter set in the T-D plane:

figure(1)
plot(T,D,'*r','MarkerSize',20)
set(gca, 'fontsize',16)

% We also plot the two axes, and the parabola of equation D = .25*T^2,
%    since the qualitative behavior of the system is determined by
%    the position of (T,D) w.r.t. these 3 lines:

hold on
min_T = -4;
max_T =  4;
min_D = -1;
max_D =  4;
plot([min_T max_T],[0 0])
plot([0 0],[min_D max_D])
T_axis = min_T:.01:max_T;
plot(T_axis,.25*T_axis.^2)
hold off
title(['a = ' num2str(a) '    b = ' num2str(b) '    c = ' num2str(c) '    d = ' num2str(d) '    e = ' num2str(e) '    f = ' num2str(f) ])
grid
xlabel('TRACE (= a + d)')
ylabel('DETERMINANT (= a*d - b*c)')

 
%% DIRECTION FIELD 

min_x = -10;
max_x =  10;
min_y = -10;
max_y =  10;
% min_x =  0;
% max_x =  16;
% min_y = -3;
% max_y =  13;

% open new figure:
figure(2)

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

% Compute the derivatives at each grid point:
dx = F(xx,yy);
dy = G(xx,yy);

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

% 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)

% find the equilibrium:

if ~(D == 0)
    
    x0 = (b*f - e*d)/D;
    y0 = (c*e - a*f)/D;
    
    % if the eigenvalues are real-valued, plot the separatrices:
    
    if isreal(eigenvalue_1)
        separatrix_1 = @(x,y) eigenvector_1(2)*(x - x0) - eigenvector_1(1)*(y - y0);
        h_sep_1 = ezplot(separatrix_1,[min_x,max_x,min_y,max_y]);
        separatrix_2 = @(x,y) eigenvector_2(2)*(x - x0) - eigenvector_2(1)*(y - y0) ;
        h_sep_2 = ezplot(separatrix_2,[min_x,max_x,min_y,max_y]);
        set(h_sep_1,'color','k')
        set(h_sep_2,'color','k')
    end
    
end


%% 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.

dt = .0001;      % time step
t_max = 50;     % length of simulation
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

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(2)
    X(:,1) = ginput(1);    % click mouse on desired initial state
    i = 1;
    while i <= n_time_steps && X(1,i)>min_x && X(1,i)<max_x && X(2,i)>min_y && X(2,i)<max_y
        i = i + 1;
        X(:,i) = X(:,i-1) + (A*X(:,i-1) + non_homogeneous_term)*dt;
    end
    ii = 1:i-1;             % time axis for current trajectory
    plot(X(1,ii),X(2,ii),'r','linewidth',2)
    
    if ~ isempty(ii)
        figure(3)
        plot(time_axis(ii),X(1,ii),time_axis(ii),X(2,ii),'linewidth',2)
        legend('x(t)','y(t)')
        xlabel('t')
        grid
        set(gca, 'fontsize',16)
    end
    
end

title(['a = ' num2str(a) '    b = ' num2str(b) '    c = ' num2str(c) '    d = ' num2str(d) '    e = ' num2str(e) '    f = ' num2str(f) ])

if (~(D == 0)) && isreal(eigenvalue_1)
    legend([h_V_nc h_H_nc h_sep_1],'V-nc','H-nc','sep.')
else
    legend([h_V_nc h_H_nc],'V-nc','H-nc')
end

hold off
