function B=expAt(A,t)
%-----------------------------------------------------------------------
% Find exp(A*t) in the correct mathematical sense.
% Note that Matlab's exp(A*t) is an element-wise operation.
% Variables:
%   A     ... nxn square matrix
%   t     ... time (scalar)
%   done  ... keep track of calculation of Jordan normal form
%   V     ... similarity transform matrix:  A*T=T*J
%   J     ... Jordan normal form
%   expJt ... nxn matrix of exp(J*t)
%   B     ... exp(A*t)
% Algorithm: similarity transform A*T=T*J     exp(A*t)*T=T*exp(J*t)
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

      global V J done

% Find the Jordan normal form J and similarity transform matrix V
      if (isempty(done) == 1 | done==0)
        [V,J] = jordan(A);
        done=1;
      end

% Diagonal elements of exp(J*t)
      for i = 1 : size(J,1)
        expJt(i,i)=exp(J(i,i)*t);
      end

% Off-diagonal elements of exp(J*t)
      for i = 1 : size(J,1)-1
        if J(i,i+1)~=0 expJt(i,i+1)=t*exp(J(i,i)*t); end
      end

% Form exp(A*t)
      B=V*expJt*V^(-1);