Number + matrix: behavior and the resulting problems

Hello,

I've noticed that if you have a matrix and add a number to it, it adds the
number to each element in the matrix.

>From my engineering courses we've always used the convention that
number+matrix is actually interpreted as number*Identity + matrix. This
would be the way to solve it, except Maxima's solve doesn't like

solve([
unknown_matrix^^0 . ident(n) = linear_eq1(B1, B2, ..., Bn),
unknown_matrix^^1 . ident(n) = linear_eq2(B1, B2, ..., Bn),
...
unknown_matrix^^(n-1) . ident(n) = linear_eqn(B1, B2, B3, ..., Bn)],
[B1, B2, ..., Bn]);

If I drop the ident(n)

solve([
unknown_matrix^^0 = linear_eq1(B1, B2, ..., Bn),
unknown_matrix^^1 = linear_eq2(B1, B2, ..., Bn),
...
unknown_matrix^^(n-1) = linear_eqn(B1, B2, B3, ..., Bn)],
[B1, B2, ..., Bn]);

Then this returns the wrong answer as unknown_matrix^^0 reduces to 1, then
it is added/subtracted from each element. Even when I use
declare(unknown_matrix, nonscalar)$
unknown_matrix^^0;
This returns the value 1, which is scalar. I'm not sure what the
appropriate behavior should be, but the current behavior feels wrong to me.

So I'm wondering if there is anyway to make maxima do this:
variable^^0 = symbolic_identity,
where the symbolic identity behaves like an identity matrix, whose size is
determined at time of numerical operation (i.e., when it is added to a
square matrix of numbers or variables)

I have ugly workarounds, but I'm looking for an elegant solution.