/*
Compute the exponential of a matrix using the Putzer Algorithm.
Copyright (C) 2003  Dan Stanger

This library is free software; you can redistribute it and/or modify it
under the terms of the GNU Library General Public License as published
by the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.

This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

Dan Stanger dan.stanger@ieee.org
*/

Putzer(a,eigenvalues,var):=block([n:length(eigenvalues),p,i,r,lambda],
   i:ident(n),
   local(r,p,lambda),
   array(r,n),
   array(p,n),
   array(lambda,n),
   /* listarray(lambda,eigenvalues), */
   for k:1 thru n do lambda[k]:eigenvalues[k],
   p[0]:i,
   p[1]:a-lambda[1]*i,
   r[1]:rhs(ic1(
      ode2('diff(arraymake(r,[1]),var)-lambda[1]*arraymake(r,[1]),
           arraymake(r,[1]),var),
        var=0,arraymake(r,[1])=1)),
   for j:2 thru n do (
      r[j]:rhs(ic1(
         ode2('diff(arraymake(r,[j]),var)-lambda[j]*arraymake(r,[j])-r[j-1],
              arraymake(r,[j]),var),
           var=0,arraymake(r,[j])=0)),
      p[j]:i,for k:1 thru j do p[j]:p[j].(a-lambda[k]*i)
   ),
   sum(r[j+1]*p[j],j,0,n-1)
)$
/* Putzer(matrix([3,-1],[1,1]),[2,2],t) */
/* Putzer(matrix([1,1,0],[0,1,0],[-1,2,2]),[1,1,2],t) */
/* Putzer(matrix([2,1,0,0],[0,2,1,0],[0,0,2,0],[0,0,0,3]),[2,2,2,3],t) */