%% Case of complex (nonreal) eigenvalues r_1=\alpha+i\beta, r_2=\alpha-i\beta

clearvars

% Example: consider the ODE system \vec y\prime= \left[\matrix{ 4 & 5 \cr -4 &0} \right]\vec y

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

% Find the eigenvalues: solve the characteristic equation
%         \det(A-rI)=0 \\
% r^2-4r+20=0
% This gives the eigenvalues
% r_1 = 2+4i,\;r_2 = 2-4i. Note that r_2 = \overline{r_1}.

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

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

% For r_1 =2-4i we get the eigenvector \vec v^{(1)}=
% \left[\matrix{
% -\frac12+i \cr 1
% }\right]. This gives the complex solution
%         e^{r_1 t}\,\vec v^{(1)}=e^{(2-4i)t} \left[\matrix{
% -\frac12+i \cr 1
% }\right]
% For r_2 =2+4i we get the eigenvector \vec v^{(2)}=\left[\matrix{
% -\frac12-i \cr 1
% }\right]. This gives the complex solution
%         e^{r_2 t}\,\vec v^{(2)}=e^{(2+4i)t} \left[\matrix{
% -\frac12-i \cr 1
% }\right]
%% RECIPE: Find a fundamental set \vec Y^{(1)} (t), \; \vec Y^{(2)}(t)  of real solutions
% Pick a complex solution \vec y(t) = e^{(\alpha+i\beta)t}\left(\vec a+i\,\vec b\right)
% Let \vec Y^{(1)}(t)=\text{Re} \vec y(t) = e^{\alpha t}\left(
% \cos(\beta t)\vec a-\sin(\beta t)\vec b
% \right)
% Let \vec Y_2(t)=\text{Im} \vec y(t) = e^{\alpha t}\left(
% \cos(\beta t)\vec b+\sin(\beta t)\vec a
% \right)
% Here we have r=2+4i and \vec v =
%  \left[\matrix{-\frac12 \cr 1}
% \right] +
% i\, \left[\matrix{-1 \cr 0}
% \right]. This gives the fundamental set of solutions
% \vec Y_1 (t) =
% e^{2t}  \left( \cos(4t)  \left[\matrix{-\frac12 \cr 1}
% \right]
% +\sin(4t)\left[\matrix{1 \cr 0}
% \right]
% \right)
% \vec Y^{(2)}(t) =
% e^{2t}  \left( \cos(4t)
%
% \left[\matrix{-1 \cr 0}
% \right]
%
%
% +\sin(4t)
%  \left[\matrix{-\frac12 \cr 1}
% \right]
%
% \right)

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
Y1 = subs(ysol,{C1,C2},{1,0})
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),\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(Y2(1),Y2(2),[-2,2],'k');              % plot trajectory for Y2 in black
fplot(-Y1(1),-Y1(2),[-2,2],'b:');           % plot trajectory for -Y1 in blue, dotted
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')