%% Solving the ODE system \vec y\prime = A\vec y
% We are given a matrix A = \left[
% \matrix{a_{11} & a_{12} \cr a_{21} & a_{22} } \right] and want to find \left[ \matrix{ y_1(t) \cr y_2(t) }
% \right] such that
%          \left[ \matrix{ y_1\prime(t) \cr y_2\prime(t) }
% \right] =
% \left[
% \matrix{a_{11} & a_{12} \cr a_{21} & a_{22} } \right]
% \left[ \matrix{ y_1(t) \cr y_2(t) }
% \right]
% Example: Solve the initial value problem
%          \left[ \matrix{ y_1\prime(t) \cr y_2\prime(t) }
% \right] =
% \left[
% \matrix{1 & 2 \cr 4 & 3} \right]
% \left[ \matrix{ y_1(t) \cr y_2(t) }
% \right]  with initial condition
% \left[ \matrix{ y_1(0) \cr y_2(0) }
% \right]
% = \left[ \matrix{ 1 \cr 0 }
% \right]
%% Key idea: use the guess \vec y(t) = e^{rt} \,\vec v
% We need to find a number r and a vector \vec v\ne \vec 0 such that \vec y\prime=A\vec y:
% We obtain the equation A\vec v = r \vec v. With the identity matrix I= \left[
% \matrix{1 & 0 \cr 0 & 1}
% \right] we can write this as
%          (A-rI)\vec v=\vec 0
% Since \vec v\ne 0 the matrix A-rI must be singular:
%        \det(A-rI)=\det \left[ \matrix{
% a_{11}-r&a_{12} \cr a_{21} & a_{22}-r
% }\right] =0
% This gives the characteristic equation
%        r^2 - \text{trace}\cdot r+ \text{determinant} = 0
% where \text{trace} = a_{11}+a_{22}=4 and \text{determinant}=\det A=-5:
%        r^2-4r-5=0
% Solving this quadratic equation gives the eigenvalues r_1=5,\; r_2=-1

A = sym([1 2; 4 3])
eig(A)            % find eigenvalues of matrix A

% For each eigenvalue r we need to find eigenvectors \vec v satisfying (A-rI)\vec v = \vec 0:
%          \left[\matrix{ a_{11}-r & a_{12} \cr a_{21} & a_{22}-r
% }\right]
% \left[\matrix{ v_1 \cr v_2
% }\right]
% =\left[\matrix{ 0 \cr 0
% }\right]
% For the eigenvalue r=5 we obtain the linear system
%          \left[\matrix{ -4 & 2 \cr 4 & -2
% }\right]
% \left[\matrix{ v_1 \cr v_2
% }\right]
% =\left[\matrix{ 0 \cr 0
% }\right]
% Note that the matrix is singular, so we can find nonzero solution vectors. We obtain -4v_1+2v_2 = 0 or 2v_1 = v_2, so
%          \vec v = \left[\matrix{ 1 \cr 2} \right]
% is an eigenvector for r=5 (multiplying this vector by any nonzero number would also be an eigenvector).
% For the eigenvalue r=-1 we obtain the linear system
%          \left[\matrix{ 2 & 2 \cr 4 & 4
% }\right]
% \left[\matrix{ v_1 \cr v_2
% }\right]
% =\left[\matrix{ 0 \cr 0
% }\right]
% Note that the matrix is singular, so we can find nonzero solution vectors. We obtain 2v_1+2v_2 = 0 or v_1 = -v_2, so
%          \vec v = \left[\matrix{ 1 \cr -1} \right]
% is an eigenvector for r=-1 (multiplying this vector by any nonzero number would also be an eigenvector).
% Find eigenvalues and eigenvectors with Matlab: The command [V,D]=eig(A) gives
% a matrix V = \left[
% \vec v^{(1)}, \vec v^{(2)} \right] containing the eigenvectors
% a diagonal matrix D = \left[ \matrix{ r_1 & 0 \cr 0 & r_2 } \right] containing the eigenvalues

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

% Matlab gives
% the eigenvalue r_1=-1 with eigenvector \vec v = \left[\matrix{-1 \cr 1}\right]
% the eigenvalue r_2=5 with eigenvector \vec v = \left[\matrix{\frac12 \cr 1}\right]
% The general solution of the ODE \vec y\prime=A\vec y is now given by
%         \vec y(t) =
% c_1 e^{r_1 t} \,\vec v ^{(1)}
% +c_2 e^{r_1 t} \,\vec v ^{(2)}
% Here we obtain
%          \vec y(t) = c_1 e^{-t} \left[ \matrix{-1 \cr 1 } \right]
% +  c_2 e^{5t} \left[ \matrix{\frac12 \cr 1 } \right]
% For the initial condition we plug in t=0 and obtain the equation
%         c_1 \left[ \matrix{-1 \cr 1 } \right]
% +  c_2  \left[ \matrix{\frac12 \cr 1 } \right] =
% \left[\matrix{ 1 \cr 0 }\right]     or    \left[\matrix{
% -1 & \frac12 \cr 1 & 1
% }\right]
% \left[\matrix{
% c_1 \cr c_2
% }\right] =
% \left[\matrix{
% 1 \cr 0
% }\right]
% We can solve this linear system in Matlab as follows

c = V\[1;0]                % solve linear system  V c = [1;0]

% Hence the solution of the initial value problem is
%          \vec y(t) = -\frac23 e^{-t} \left[ \matrix{-1 \cr 1 } \right]
% +  \frac23 e^{5t} \left[ \matrix{\frac12 \cr 1 } \right]
%% Solving an ODE system with dsolve
% We have to tell Matlab that y(t) is an array of size [2,1]:

syms y(t) [2 1]             % define vector valued function y(t) = [y1(t);y2(t)]
syms C1 C2

%% Find general solution C_1 \vec Y^{(1)}(t)+C_2 \vec Y^{(2)}(t)

sol = dsolve( diff(y)==[1 2;4 3]*y);
ysol = 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[-3,3]

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

%
%% Solve initial value problem
% We have the initial condition \vec y(0)= \left[\matrix{1 \cr 0} \right]:

sol = dsolve( diff(y)==[1 2;4 3]*y, y(0)==[1;0] )   % solve initial value problem
ysol = subs(y(t),sol)       % in y(t) substitute the values from sol

%% Plot vector field and trajectory in the phase plane
% You need to download the file vectorfield.m and put it in the directory where your mlx files are.
% We plot the curve \vec y(t) in the y_1,y_2 plane for t\in[-3,3]. The inital point \vec y^{(0)}=\left[\matrix{1\cr 0}\right] for t=0 is marked with a green dot.

f = @(t,y) A*y;                            % define function f(t,y) for vectorfield
vectorfield(f,-8:8,-8:8); hold on
plot(1,0,'g.','MarkerSize',20)             % mark initial point as green dot
fplot(ysol(1),ysol(2),[-3,3],'g');         % plot trajectory for ysol as green curve
hold off; axis equal; axis([-8 8 -8 8]); 
xlabel('y_1'); ylabel('y_2')
title('solution of initial value problem (green curve)')

%%
%
%