% script hw13_20
% calculates the critical points of a phase plane system,
% evaluates the Jacobian, and computes the eigenvalues and 
% eigenvectors at each critical point
function hw13()

clear all
% find critical points 
% # 20
% [x0,y0] = solve('6*x0 - 5*y0 +x0^2 = 0','2*x0-y0+y0^2 = 0');
% xp = sym('6*x - 5*y +x^2'); yp = sym('2*x-y+y^2');
% #21
% [x0,y0] = solve('x0 + 2*y0 + x0^2 + y0^2 = 0','2*x0 - 2*y0 + 3*x0*y0 = 0');
% xp = sym('x + 2*y + x^2 + y^2'); yp = sym('2*x - 2*y + 3*x*y');
% #30
%[x0,y0] = solve('y0 - 1 = 0',' - y0 + x0^2 = 0');
%xp = sym('y - 1'); yp = sym('- y + x^2');
% #32
 [x0,y0] = solve('x0*y0 - 2 = 0',' x0 - 2*y0 = 0');
 xp = sym('x*y - 2'); yp = sym('x - 2*y');
simplify(x0)
simplify(y0)
% evaluate symbolic expressions to produce numerical values
x1 = eval(x0);
y1 = eval(y0);
% define symbolic functions for F and G
% compute and print Jacobian
Jacobian = jacobian([xp;yp])
n = length(x0)
for i = 1:n
    sprintf('------------------------')
    i
% set x1, y1, the symbolic variables used in the Jacobian,
% to one of the critical points
    x = x1(i);y = y1(i);
    [x,y]
% evaluate the components of the Jacobian
    a = eval(Jacobian(1,1));
    b = eval(Jacobian(1,2));
    c = eval(Jacobian(2,1));
    d = eval(Jacobian(2,2));
% define the linearization matrix at i-th critical point
% and compute its eigenvectors V and eigenvalues D
    A = [a, b; c, d]
    [V,D]=sys_symb(A)
end


%   try to compute eigenvectors with integer entries

function [V,D] = sys_symb(A)
[m,n] = size(A);
err = .5*ones(m,1);
if m ~= n 
   fprintf('matrix must be square')
   break
end
[V,D] = eig(A); V1 = V;
% scale eigenvectors to reveal integral coefficient forms (if possible)
for i = 1:m
   vmax = abs(V(1,i));
   for j = 2:m
      vmax = max(vmax,abs(V(j,i)));
   end    
   if vmax == 0
      fprintf('defective eigenproblem (sorry!)')
   end
%  scale i-th eigenvector by its absolute-maximum component
   V(:,i) =  V(:,i)/vmax;
   kk = 0; j=1;
   while kk == 0 & j <= 10000
% multiply the i-th eigenvector by j (which can become as large as 10000 here)
% and check how close to an integer are the resulting components
      V1(:,i) = V(:,i)*j;
%  this adds up the distance of each component to the nearest integer;
%  if it is acceptably small, then this scaled eigenvector is accepted 
%  as the  desired integer-component form sought
      if sum(abs(abs(V1(:,i)) - floor(abs(V1(:,i))+err))) <= .0001
         V(:,i) = V1(:,i);
         kk = 1;
      end
      j = j+1;
   end
end
