%% Solution Exam 3 (MATH246, Spring 2024)
%% Problem 1(a),(b)

syms x y
f = 1-2*y;
g = x^2-y^2;

% Find stationary points:

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

% We found two stationary points: (-\frac12,\frac12) and (\frac12,\frac12).

J = jacobian([f;g],[x;y])
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

% Result:
% (-\frac12,\frac12) is saddle (unstable, not repelling)
% (\frac12,\frac12) is spiral sink (stable, attracting)
%% Problem 1(c)
% The stationary points for the nonlinear system have the same type and stability as in 1(b) (only for the center we may have differences for the nonlineary ODE system).

f = @(t,y) [ 1-2*y(2) ; y(1)^2-y(2)^2 ];     % we now use y1,y2 instead of x,y
for a1 = -1:.25:1
  for a2 = 0:.125: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,-1:.25:1,0:.125:1); 
hold off; axis equal; axis([-1 1 0 1]);

%
%
%% Problem 2

A = sym([6 2; 5 9])
[eigenvectors,eigenvalues] = eig(A)

% We have two different positive eigenvalues, so this is a nodal source (unstable, repelling).

phaseportrait(A)

%% Problem 3

A = sym([0 -2; 8 0])
eigenvalues = eig(A)   

% We have two imaginary eigenvalues, so this is a center (stable, not attracting), counterclockwise.

phaseportrait(A)

%
%% 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