%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Math 3520: Differential Equations and Dynamical Systems
%   Linear first-order systems: X'=AX, 2x2 matrix A
%   Flow fields for various classes of eigensystems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% DO THIS FIRST:
% Initialize the grid in phase space where arrows will
% indicate the flow X'=(x',y') at position X=(x,y):
t=-4:4; [x,y]=meshgrid(t,t);  % 2-D grid 
dt=0.1; % Small increment for short arrows

% NOTE: The plotting function is the same for all examples.
% dx = A(1,1)*x*dt + A(1,2)*y*dt;
% dy = A(2,1)*x*dt + A(2,2)*y*dt;
% quiver(x,y,dx,dy);

%%% DISTINCT REAL EIGENVALUES:
% Both negative:  lambda1 < lambda2 < 0
A=[-1 0; 0 -2]; A
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);

% Both positive:  0 < lambda1 < lambda2
A=[ 1 0; 0 2 ]; A
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);

% One negative, one positive:  lambda1 < 0 < lambda2
A=[ -1 0; 0 2 ]; A
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);

%%% DISTINCT COMPLEX EIGENVALUES:
% Purely imaginary: Re(lambda) = 0
A=[ 0 2; -2 0 ]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);

% Positive real part: Re(lambda) > 0
A=[ 1 2; -2 1 ]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);

% Negative real part: Re(lambda) < 0
A=[ -1 2; -2 -1 ]; A, eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);


%%% REPEATED REAL EIGENVALUES, DEGENERATE EIGENVECTORS:
% Positive, degenerate: lambda1 = lambda2  > 0
A=[ 1 1; 0 1 ]; A, [eVec,eVal]=eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);

% Negative, degenerate: lambda1 = lambda2  < 0
A=[ -1 1; 0 -1 ]; A, [eVec,eVal]=eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);

% Zero, degenerate: lambda1 = lambda2  = 0
A=[ 0 1; 0 0 ]; A, [eVec,eVal]=eig(A)
dx = A(1,1)*x*dt + A(1,2)*y*dt;
dy = A(2,1)*x*dt + A(2,2)*y*dt;
quiver(x,y,dx,dy);