Linear Systems



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?