One missing maxima feature, that I have the math behind, but not a program

Hello all,

I was recently making use of the matrixexp function, when I realized there
is a related operation that Maxima does not yet appear to support natively.

matrixexp(A * t) is often useful for solving linear differential equations.
Similarly, A^^t is useful for solving (discrete) difference equations.

However, I could not find an operation to perform A^^t. So anyways I threw
together a sample maxima program which carries out a matrix exponential and
an A^^t by "hand" (it isn't automated). I was wondering if anyone cared to
convert the A^^t to some method that Maxima can call?

The maxima file is included as text below (I can send the file to anyone
who wants it):
---------------------------------------------
/* [wxMaxima: input   start ] */
m1: matrix([1, 1, 0, 0],[0, 1, 1, 0],[0,0,1,-1/8],[0,0,1/2,1/2]);
matrixexp(m1*t);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eivals(m1);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eq1: exp(A*t) = B11 * exp(3/4 * t) + B12 * t * exp(3/4 * t) + B21 * exp(t)
+ B22 * t * exp(t);
eq2: facsum( diff(eq1,t,1), B11, B12, B21, B22);
eq3: facsum( diff(eq1,t,2), B11, B12, B21, B22);
eq4: facsum( diff(eq1,t,3), B11, B12, B21, B22);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
declare(I,nonscalar)$
declare(A,nonscalar)$
declare(B11,nonscalar)$
declare(B12,nonscalar)$
declare(B21,nonscalar)$
declare(B22,nonscalar)$
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eq5: solve( [lhs(eq1) . I = rhs(eq1), eq2, eq3, eq4], [B11,B12,B21,B22] ),
t=0;
slns: eq5, A^3=A^^3, A^2=A^^2, A = m1, I = ident(4);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eq6: radcan(eq1), slns[1][1], slns[1][2], slns[1][3], slns[1][4];
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
is( matrixexp(m1*t) = rhs(eq6));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
/*                       */
/* Now the discrete case */
/*                       */
m1;
m1^^2;
m1^^3;
m1^^4;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
m1: matrix([1, 1, 0, 0],[0, 1, 1, 0],[0,0,1,-1/8],[0,0,1/2,1/2]);
declare(I,nonscalar)$
declare(A,nonscalar)$
declare(B11,nonscalar)$
declare(B12,nonscalar)$
declare(B21,nonscalar)$
declare(B22,nonscalar)$
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eq1: A^^t = B11 * (3/4) ^ t + B12 * t * (3/4)^t + B21 * 1^t + B22 * t *
1^(t);
eq2: eq1, t=w+1$
eq2: eq2, w=t;
eq3: eq2, t=w+1$
eq3: eq3, w=t;
eq4: eq3, t=w+1$
eq4: eq4, w=t;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
expand(-(48*I-64*A^^3+176*A^^2-160*A)/3);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eq5: solve( [lhs(eq1) . I = rhs(eq1), eq2, eq3, eq4], [B11,B12,B21,B22] ),
t=0;
slns: eq5, A^3=A^^3, A^2=A^^2, A = m1, I = ident(4);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eq6: radcan(eq1), slns[1][1], slns[1][2], slns[1][3], slns[1][4];
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
eq6, t=2;
/* [wxMaxima: input   end   ] */