%% Phase portraits for autonomous linear ODE system
% We consider the ODE system
%         \vec y\prime(t) = A\vec y(t)
% with a nonsingular matrix A\in \mathbb{R}^{2\times2}. In this case there is one stationary point at \vec 0.
% We obtain the following types of stationary points:
%% REAL EIGENVALUES:
% both positive
% different
% nodal source
% unstable
%
% equal
% radial source or twist source (deficient case)
% unstable
% both negative
% different
% nodal sink
% stable
%
% equal
% radial sink or twist sink (deficient case)
% stable
% positive & negative
%
% saddle
% unstable
% Equal eigenvalues: if there are two eigenvectors we get a radial sink/source. If there is only one eigenvector we get a twist sink/source.
%% NONREAL EIGENVALUES:
% real part positive
% spiral source
% unstable
% real part negative
% spiral sink
% stable
% real part zero
% center
% stable
% We get a general solution \vec y(t) = C_1 \vec y ^{(1)}(t)+C_2 \vec y^{(2)}(t)
% In the phase portaits below we always show eight trajectories: We use C_1 =-1,0,1 and C_2 =-1,0,1 (for C_1 =C_2=0 we just get the constant solution \vec y(t)=\vec 0).

clearvars

%% Two different negative eigenvalues: NODAL SINK

A = sym([-4 -3; -1 -2])

% The eigenvalues are r_1=-5,\;r_2=-1

[V,D] = eig(A)
phaseportrait(A)

%% Note: for two different POSITIVE eigenvalues we get a NODAL SOURCE (phase portrait as above, but with arrows pointing in the other direction).
%% positive and negative eigenvalue: SADDLE

A = sym([1 2;4 -1])

% The eigenvalues are r_1=3,\;r_2=-3

[V,D] = eig(A)
phaseportrait(A)

%% two equal negative eigenvalues, two eigenvectors: RADIAL SINK

A = sym([-2 0;0 -2])

% The eigenvalues are r_1=r_2=-2

[V,D] = eig(A)
phaseportrait(A)

%% Note: for two equal POSITIVE eigenvalues, two eigenvectors we get a RADIAL SOURCE (phase portrait as above, but with arrows pointing in the other direction).
%% two equal negative eigenvalues, one eigenvector: TWIST SINK

A = sym([-1 -1;1 -3])

% The eigenvalues are r_1=r_2=-2. There is only one eigenvector \vec v = \left[
% \matrix{1\cr 1}
% \right]:

[V,D] = eig(A)

% We have to find a generalized eigenvector \vec w satisfying (A-rI)\vec w  = \vec v: We can use e.g. \vec w = \left[
% \matrix{1\cr 0}
% \right]:

[V,D] = jordan(A)   % find eigenvectors and generalized eigenvectors
phaseportrait(A)

%% Note: for two equal POSITIVE eigenvalues, only one eigenvector we get a TWIST SOURCE (phase portrait as above, but with arrows pointing in the other direction).
%% nonreal eigenvalues, real part negative: SPIRAL SINK

A = sym([-1 5;-2 -3])

% The eigenvalues are -2\pm 3i.

[V,D] = eig(A)
phaseportrait(A)

%% Note: for nonreal eigenvalues with POSITIVE real part we get a SPIRAL SOURCE (phase portrait as above, but with arrows pointing in the other direction).
%% nonreal eigenvalues, real part zero: CENTER

A = sym([1 -5; 1 -1])

% The eigenvalues are \pm 2i.

[V,D] = eig(A)
phaseportrait(A)

%
%% Function for drawing phase portrait

function phaseportrait(A)
syms y(t) [2 1]
syms C1 C2
sol = dsolve(diff(y)==A*y);             % solve ODE 
ysol = subs(y(t),sol);
f = @(t,y) A*y;                         % define function f(t,y) for vectorfield
vectorfield(f,-5:5,-5:5); hold on
for c1=[-1 0 1]
  for c2=[-1 0 1]
    ys = subs(ysol,{C1,C2},{c1,c2});    % use C1=c1, C2=c2
    fplot(ys(1),ys(2),[-5 5]);          % plot trajectory for t=-5...5
  end
end
hold off; axis equal; axis([-5 5 -5 5])
end