%-----------------------------------------------------------------------
% Find the Jordan normal (canonical) form
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

      global V J done

% Program Header -------------------------------------------------------
      disp('Find the Jordan normal form of the following matrix')
      disp('     [ 3 -1.5 3]')
      disp('  A= [ 0  2   2]')
      disp('     [ 0 -0.5 4]')
      disp(' ')

% Provide matrix A------------------------------------------------------
      A= [ 3 -1.5 3
           0  2   2
           0 -0.5 4 ];

% Output results--------------------------------------------------------
      expat(A,0);
      disp('Jordan normal form is:')
      disp(J)
      disp('Similarity transform matrix is:')
      disp(V)
      disp('Check: V*J*V^(-1):')
      disp(V*J*V^(-1))
      disp('It should be identical to A:')
      disp(A)

%-----------------------------------------------------------------------
% Example of the use of exp(A*t) in solvind dx/dt=exp(A*t)*x0
% Initial condition
      x0=[1; 2; 0.5];
% Integrate dxdt=A*x
      for i=1:100
        t(i)=0.01*i;
        x(:,i)=expat(A,t(i))*x0;
      end

% Plots-----------------------------------------------------------------
% plot x vs. t
      subplot(2,2,1)
%     plot(t,x(1,:), t,x(2,:), t,x(3,:))    %OK
      plot(t,x)
      xlabel('Time'), ylabel('States'), title('Dynamic States')
      legend('x_1', 'x_2', 'x_3', 2)

% plot x1 vs. x2
      subplot(2,2,2)
      plot(x(2,:),x(1,:))
      xlabel('x_2'), ylabel('x_1'), title('Phase Diagram')

% plot x1 vs. x3
      subplot(2,2,3)
      plot(x(3,:),x(1,:))
      xlabel('x_3'), ylabel('x_1'), title('Phase Diagram')

% plot x2 vs. x3
      subplot(2,2,4)
      plot(x(3,:),x(2,:))
      xlabel('x_3'), ylabel('x_2'), title('Phase Diagram')