%% Using Matlab for Autonomous Systems

clearvars; clf

% You need to download the file vectorfield.m and put it in the directory where your mlx files are.
% Example: Consider the following autonomous system:
%             \left[ \matrix{y_1\prime \cr y_2\prime}
% \right]
% =
% \left[\matrix{y_1+y_1 \,y_2 \cr
% y_2-y_1\,y_2-1}
% \right]  \\
% \vec y\prime = \vec f(\vec y)
%% Define the right hand side function f

f = @(t,y) [ y(1)+y(1)*y(2) ; y(2)-y(1)*y(2)-1 ];

% Note that you must use @(t,y) ..., despite the fact that the functions do not depend on t.
%% Plot the vector field
% for y_1 between -1 and 3 and y_2 between -2 and 2 using steps of size .5:

vectorfield(f,-1:.5:3,-2:.5:2); xlabel('y_1'); ylabel('y_2')

%
%% Plot trajectories for many initial points \left[ \matrix{a_1 \cr a_2}
% \right]
% We can solve the ODE with initial condition \vec y(0)=\vec a for t going from 0 to 4 and plot the trajectory using
%        [ts,ys] = ode45(f,[0,4],[a1;a2]);
%        plot(ys(:,1),ys(:,2))
% We can solve this initial value problem t going from 0 to 4 and plot the trajectory using
%        [ts,ys] = ode45(f,[0,-4],[a1;a2]);
%        plot(ys(:,1),ys(:,2))
% The ode45 commands may generate a number of warning messages since the solutions does not exist for t in the whole interval [-4,4]). We can turn off these warnings with the command

warning('off','MATLAB:ode45:IntegrationTolNotMet')

% Now we plot the trajectories going through the points \left[ \matrix{a_1 \cr a_2}
% \right] with a_1=-1,-.5,\ldots,3 and y_2=-2,-1.5,\ldots,2

for a1=-1:.5:3
  for a2=-2:.5:2
    [ts,ys] = ode45(f,[0,4],[a1;a2]); 
    plot(ys(:,1),ys(:,2),'c'); hold on
    [ts,ys] = ode45(f,[0,-4],[a1;a2]);
    plot(ys(:,1),ys(:,2),'c')
  end
end
vectorfield(f,-1:.5:3,-2:.5:2);         % plot vector field
hold off; axis([-1 3 -2 2]); xlabel('y_1'); ylabel('y_2')

%
%% Investigation of stationary points
% Define symbolic variables

syms t y1 y2 

% Define symbolic expression F by using f with the vector [y1;y2]:

F = f(t,[y1;y2])

%% Find the stationary points
% Find \vec y where \vec F = \vec 0 using solve:

[y1s,y2s] = solve(F,[y1;y2])   % find all solutions (y1,y2) of F=[0;0]
m = length(y1s)                 % number of solutions found

% We have found two stationary points: (0,1) and (2,-1)
%% Find the Jacobian matrix:

J = jacobian(F,[y1;y2])

%% Investigate the stationary points
% Result:
% stationary point (0,1): eigenvalues 1,2, nodal source
% stationary point (2,-1): eigenvalues -2,1, saddle

for k=1:m
  stationary_point = [y1s(k),y2s(k)]

% Evaluate the Jacobian at the stationary point:

  A = subs(J,{y1,y2},{y1s(k),y2s(k)})
  [eigvect,eigvals] = eig(A)
  figure
  phaseportrait(A); 
  title(sprintf('linearized problem for stationary point (%g,%g)',stationary_point))
end

%
%% 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]); xlabel('y_1'); ylabel('y_2')
end