Given a coefficient matrix and a right-hand-side vector, this function
will solve for x_1, ..., x_n:
msolve(A, b) := block([n, vars, x, eqns],
n: length(transpose(A)),
vars: makelist(concat(x, i), i, 1, n),
x: covect(append(vars, [1])),
eqns: flatten(args(echelon(addcol(A, b)) . x)),
solve(eqns, vars));
For instance:
(%i18) A: matrix([1,1,-2], [3,-1,-7], [1,1,1], [2,2,-4])$
(%i19) b: covect([-3,2,0,-6])$
(%i20) msolve(A, b);
(%o20) [[x1 = - 2, x2 = 3, x3 = - 1]]
Problem is, it seems a little ad-hoc; does Maxima have a more elegant,
native facility?