%% Solution of Practice Problem for Exam 3

clearvars

% Consider the ODE system
%          \left[ \matrix{
%        x\prime \cr  y\prime
% }\right]
% =
% \left[ \matrix{
%    4x+y+xy \cr x+4y+y^2
% }\right]

syms x y
f = 4*x+y+x*y;
g = x+4*y+y^2;

%% Find all stationary points:
% Find \left[ \matrix{
%        x \cr  y
% }\right]  with \left[ \matrix{
%        f(x,y) \cr  g(x,y)
% }\right]
% =
% \left[ \matrix{
%        0 \cr  0
% }\right] :

[xs,ys] = solve([f;g]==[0;0],[x;y])
m = length(xs)          % number of solutions we found

% We found three stationary points (0,0), (3,-3), (-5,-5).
%% Find Jacobian matrix J:

J = jacobian([f;g],[x;y])

%% For each stationary point (x_*,y_*): find matrix A=J(x_*,y_*)
%% Find eigenvalues and eigenvectors
%% For LINEARIZED PROBLEM: find type and stability of stationary point, sketch phase portrait
% Result:
% (0,0) is nodal source (unstable, repelling)
% (3,-3) is saddle point (unstable, not repelling)
% (-5,-5) is nodal sink (stable, attracting)

for k=1:m
  stationary_point = [xs(k),ys(k)]
  A = subs(J,{x,y},{xs(k),ys(k)})       % evaluate Jacobian at stationary point
  [eigenvectors,eigenvalues] = eig(A)  
  figure
  phaseportrait(A); 
  title(sprintf('linearized problem for stationary point (%g,%g)',stationary_point))
end

%
%% For nonlinear problem: What can you say about type and stability of the stationary points?
%% What does the phase portrait for the nonlinear problem look like?
% For the nonlinear problem we will have the same types and stabilities as for the linearized problem (only in the case of a center things can be different for the nonlinear problem).
% We will have the same tangent directions of the trajectories at the stationary point.
% We use y_1, y_2 in place of x,y and define the function \vec f(t,\vec y) as follows:

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

% For many initial points \vec a: use ode45 and
% solve ODE with initial condition \vec y(0)=\vec a for t going from 0 to 5, plot trajectory in phase plane
% solve ODE with initial condition \vec y(0)=\vec a for t going from 0 to -5, plot trajectory in phase plane

for a1 = -6:5
  for a2 = -6:1
    [ts,ys] = ode45(f,[0,5],[a1;a2]);  plot(ys(:,1),ys(:,2),'c'); hold on
    [ts,ys] = ode45(f,[0,-5],[a1;a2]); plot(ys(:,1),ys(:,2),'c')
  end
end
vectorfield(f,-6:5,-6:1); 
axis([-6 5 -6 1]); hold off

% We can see
% nodal source at (0,0)
% saddle at (3,-3)
% nodal sink at (-5,-5)
% and the phase portrait close to each stationary points corresponds to the phase portrait of the linearized problem.
%
%% 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