% script syss.m
% it can find the eigenvalues, eigenvectors and coefficients of the
% initial value problem for the system:
%        dx/dt  = A x ; x(0) = x0
% where x is a nX1 vector, and A is a nXn matrix
% the initial conditions can be for several IVP at once, i.e. x0 can
% be a nXm matrix, whose columns are IC for each of m different problems.
% It attempts to "rationalize the eigenvectors; if the original system
% has integer coefficients and the eigenvalues are also integers, then
% the eigenvectors can also be found with integer coefficients. This 
% script will succeed usually, and I take no responsibility if it fails!
% (try increasing the limit of 10000 trials)
%
%-----------------------------------------------------------
%  TO RUN:
%  (1) Define the matrix A
%  (2) Define the initial vector x0
%  (3) On the matlab command window type:     syss
%-----------------------------------------------------------
%
clear all
% Follow examples of various text problems of sizes 2 up to 6.
% To run, simply uncoment data of desired dimension.
%
% initialize the matrix of coefficients, A
%A = [  9 13   0   0  0 -13;...
%    -14 19 -10 -20 10   4;...
%    -30 12  -7 -30 12  18;...
%    -12 10 -10  -9 10   2;...
%      6  9   0   6  5 -15;...
%    -14 23 -10 -20 10 0];     %5.2.50
A = [139  -14  -52  -14  28;...
     -22    5    7    8  -7;...
     370  -38 -139  -38  76;...
     152  -16  -59  -13  35;...
      95  -10  -38   -7  23]; %5.4.49
%A = [13 -42 106 139;2 -16 52 70;1 6 -20 -31;-1 -6 22 33]; %5.2.46
%A = [-40 -12 54;35 13 -46;-25 -7 34]; %5.2.42
%A = [1 2; 2 1]; %5.2.1
[m,n] = size(A);
err = .5*ones(m,1);
if m ~= n 
   fprintf('matrix must be square')
   break
end
% Initialize the Initial Conditions vector, if given; else set to zero
%x0 = [1;0;-1;1;0;1]; % sample IC for 6th order system
 x0 = [1;0;-1;3;1];   % sample IC for 5th order system
%x0 = [3;0;1;3];      % sample IC for 4th order system
%x0 = [1;0;1];        % sample IC for 3rd order system
%x0 = [1;2];          % sample IC for 2nd order system
[m1,n1] = size(x0);
if m1 ~= n
   fprintf('number of initial conditions must match order of system')
   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
D
V       
C = V\x0
