From what I see maxima can only work with matrices of fixed dimensions. I cannot really believe this! I tried to compose a matrix of child matrices and invert the parent matrix and always maxima suggests to divide by the child matrices as if they are scalars. It seems declare() is only applicable for scalar types. Is this really the case? Using the matrix function works only with integer constants for the dimensions.
J:matrix([Jm, Jms],[Jsm, Js]);
x:matrix([xm], [xs]);
f:matrix([fm],[fs]);
J.x - f=0;
invert(J).f;
I find it's possible to get the results you want by declaring the variables to represent submatrices as nonscalars, and then applying a function to invert J in a way which takes noncommutative multiplication into account (invert doesn't; this is a deficiency).
Here's what I got to work.
declare ([Jm, Jms, Jsm, Js, xm, xs, fm, fs], nonscalar);
matrix_element_mult: ".";
Here is the code from the problem statement,
J:matrix([Jm, Jms],[Jsm, Js]);
x:matrix([xm], [xs]);
f:matrix([fm],[fs]);
J.x - f=0;
At this point, I'll call invert_by_lu and tell it to work in the ring of noncommutative multiplication.
invert_by_lu (J, noncommutingring);
J_inverse: expand (%);
J_inverse.f;
At this point, you could assign specific matrices to Jm, Js, Jsm, Jms, etc. and reevaluate J_inverse, via ev(J_inverse) or other means.
Note that I called expand after calling invert_by_lu. That's because a literal 1 appears in the result from invert_by_lu, which is not quite right; it should be some representation of an identity matrix. Fixing that would probably require some careful thought and some lines of code; it's not out of the question, however, I've taken the easy way out and just called expand.