%% Double eigenvalue r_1=r_2

clearvars

% There are two cases:
% We can find two linearly independent eigenvectors \vec v^{(1)}, \vec v^{(2)}
% We can find only one eigenvector \vec v(deficient case). We need to find a "generalized eigenvector" \vec w with (A-rI)\vec w = \vec v. This gives the solution \vec Y^{(2)}(t)=te^{rt}\vec v+e^{rt} \vec w.
%% Case 1: there are two linearly independent eigenvectors
% Example: consider the ODE system \vec y\prime= \left[\matrix{ 0 & 2 \cr 2 &0} \right]\vec y

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

% Find the eigenvalues: solve the characteristic equation
%         \det(A-rI)=0 \\
% (r-2)^2 = 0
% This gives the eigenvalues
% r_1 =r_2=2.

r = eig(A)      % find eigenvalues of A

% Find the eigenvectors: For each eigenvalue r find \vec v\ne \vec0 with (A-rI)\vec v = \vec 0
% Here A-rI = \left[ \matrix{ 0 & 0 \cr 0 & 0}\right] and we can find two linearly independent eigenvectors, e.g., \vec v ^{(1)}=\left[\matrix{1 \cr 0}\right], \;
% \vec v ^{(2)}=\left[\matrix{0 \cr 1}\right] (actually, any two linearly independent vectors work).

[V,D] = eig(A)  % find matrix V with eigenvectors, diagonal matrix D with eigenvalues

% Therefore we get a fundamental set of solutions
%        \vec Y^{(1)}(t) = e^{r t} \vec v^{(1)} = e^{2t} \left[\matrix{1 \cr 0}\right],
% \qquad
% \vec Y^{(2)}(t) = e^{r t} \vec v^{(2)} = e^{2t} \left[\matrix{0 \cr 1}\right]
% Solve the ODE system with dsolve:

syms y(t) [2 1]                     % define vector valued function y(t) = [y1(t);y2(t)]
syms C1 C2
sol = dsolve( diff(y)==A*y );       % find general solution
ysol = simplify( subs(y(t),sol)  )  % in y(t) substitute the values from sol

% This is the general solution, with C_1,C_2\in\mathbb{R}. We obtain \vec Y^{(1)}(t) by substituting C_1=1,\;C_2=0:

Y1 = subs(ysol,{C1,C2},{1,0})       

% We obtain \vec Y^{(2)}(t) by substituting C_1=0,\;C_2=1:

Y2 = subs(ysol,{C1,C2},{0,1})

%% Plot vector field and trajectories in the phase plane
% You need to download the file vectorfield.m and put it in the directory where your mlx files are.
% Plot the curves \pm\vec Y^{(1)}(t) and \pm\vec Y^{(2)}(t) for t\in[-2,2]

f = @(t,y) A*y;                              % define function f(t,y) for vectorfield
vectorfield(f,-20:2:20,-20:2:20); hold on
fplot(Y1(1),Y1(2),[-2,2],'b');               % plot trajectory for Y1 in blue
fplot(-Y1(1),-Y1(2),[-2,2],'b:');            % plot trajectory for -Y1 in blue, dotted
fplot(Y2(1),Y2(2),[-2,2],'k');               % plot trajectory for Y2 in black
fplot(-Y2(1),-Y2(2),[-2,2],'k:'); hold off;  % plot trajectory for -Y2 in black, dotted
axis equal; axis([-20 20 -20 20]); xlabel('y_1'); ylabel('y_2')
title('four trajectories in the phase plane')

% The general solution is \vec y(t)= e^{2t}\left[\matrix{C_1 \cr C_2}\right], so every trajectory is a straight line in the phase plane.
%
%% Case 2: there is only one linearly independent eigenvector (DEFICIENT CASE)
% Example: consider the ODE system \vec y\prime= \left[\matrix{ 5 & -4 \cr 1 &1} \right]\vec y

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

% Find the eigenvalues: solve the characteristic equation
%         \det(A-rI)=0 \\
% r^2-6r+9 = 0
% This gives the eigenvalues
% r_1 =r_2=3.

r = eig(A)      % find eigenvalues of A

% Find the eigenvectors: For each eigenvalue r find \vec v\ne \vec0 with (A-rI)\vec v = \vec 0
% Here A-rI = \left[ \matrix{ 2 & -4 \cr 1 & -2}\right] and any solution of (A-rI)\vec v=\vec0 has the form \vec v = c\left[\matrix{2 \cr 1}\right] with c\in\mathbb{R}. There is only one linearly independent eigenvector (deficient case).

[V,D] = eig(A)  % find matrix V with eigenvectors, diagonal matrix D with eigenvalues

% Note that Matlab returns only one eigenvector in V.
% This gives one solution
%        \vec Y^{(1)}(t) = e^{r t} \vec v= e^{3t} \left[\matrix{1 \cr 0}\right]
% For a second linearly independent solution we need a "generalized eigenvector" \vec w with (A-rI)\vec w = \vec v.
% We can find all eigenvectors and generalized eigenvectors with jordan(A) (find "Jordan normal form" of matrix A)

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

% The matrix D is no longer diagonal. It contains "1" above the diagonal if the corresponding vector in V is not an eigenvector, but a generalized eigenvector. Here we get
% an eigenvector \vec v = \left[\matrix{2 \cr 1}\right] satisfying A\vec v = r\vec v, yielding the solution \vec Y^{(1)}(t) = e^{r t} \vec v = e^{3t} \left[\matrix{2 \cr 1}\right]
% a generalized  eigenvector \vec w = \left[\matrix{1 \cr 0}\right] satisfying A\vec w = r\vec w + \vec v, yielding the solution \vec Y^{(2)}(t) = te^{r t} \vec v + e^{rt} \vec w
% = te^{3t} \left[\matrix{2 \cr 1}\right]
% +e^{3t} \left[\matrix{1 \cr 0}\right]
% Solve the ODE system with dsolve:

syms y(t) [2 1]                     % define vector valued function y(t) = [y1(t);y2(t)]
syms C1 C2
sol = dsolve( diff(y)==A*y );       % find general solution
ysol = simplify( subs(y(t),sol)  )  % in y(t) substitute the values from sol

% This is the general solution, with C_1,C_2\in\mathbb{R}. We obtain \vec Y^{(1)}(t) by substituting C_1=1,\;C_2=0:

Y1 = subs(ysol,{C1,C2},{1,0})       

% We obtain \vec Y^{(2)}(t) by substituting C_1=0,\;C_2=1:

Y2 = subs(ysol,{C1,C2},{0,1})

%% Plot vector field and trajectories in the phase plane
% Plot the curves \pm\vec Y^{(1)}(t) and \pm\vec Y^{(2)}(t) for t\in[-2,2]

f = @(t,y) A*y;                              % define function f(t,y) for vectorfield
vectorfield(f,-1:.1:1,-1:.1:1); hold on
fplot(Y1(1),Y1(2),[-2,2],'b');               % plot trajectory for Y1 in blue
fplot(-Y1(1),-Y1(2),[-2,2],'b:');            % plot trajectory for -Y1 in blue, dotted
fplot(Y2(1),Y2(2),[-2,2],'k');               % plot trajectory for Y2 in black
fplot(-Y2(1),-Y2(2),[-2,2],'k:'); hold off;  % plot trajectory for -Y2 in black, dotted
axis equal; axis([-1 1 -.5 .5]); xlabel('y_1'); ylabel('y_2')
title('four trajectories in the phase plane')