% Linear_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 2008

% PARAMETERS

a = 1;
b = -3;
c = 2;
d = -3;
e = 0;
f = 0;

% Matrix form:
A = [a b
     c d];
non_homogeneous_term = [e;f];
 
% THEORETICAL ANALYSIS 
 
% Eigenvalue Analysis:
[Eigenvectors,Eigenvalues] = eig(A);
eigenvalue_1 = Eigenvalues(1,1);
eigenvalue_2 = Eigenvalues(2,2);
eigenvector_1 = Eigenvectors(:,1);
eigenvector_2 = Eigenvectors(:,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)

% 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 = -6;
max_T =  6;
min_D = -5;
max_D =  5;
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

xlabel('TRACE (a + d)')
ylabel('DETERMINANT (a*d - b*c)')
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')])
    
grid

 
% DIRECTION FIELD 

% Display setup (will work in most cases, but you may have to adjust it):
min_x = -4;
max_x = 4;
step_x = .5;
min_y = -4;
max_y = 4;
step_y = .5;

% open new figure:
figure(2)

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

% Compute the derivatives at each grid point:
dx = a*x + b*y + e;
dy = c*x + d*y + f;

% 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])
axis image
xlabel('x')
ylabel('y')
hold on
grid

% 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 = .001;      % time step
t_max = 10;     % length of simulation

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

X = zeros(2,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

    X(:,1) = ginput(1);    % click mouse on desired initial weight vector

    for i = 2:n_time_steps
        X(:,i) = X(:,i-1) + (A*X(:,i-1) + non_homogeneous_term)*dt;
    end
    plot(X(1,:),X(2,:),'r')

end

hold off