*Upadate: I’ve implemented a Maxima matrix-vector equation solver with simpler Maxima-specific algorithm in a later post. That method is based on the built-in function **linsolve()**. In that post I show an example that breaks **linsolve()** but that is handled correctly by the **backsolve()** method.

Is there really not a solver in Maxima that takes matrix A and vector b and returns the solution of ? Of course we could do **invert(A).b**, but that ignores consistent systems where isn’t invertible…or even isn’t square.

Here’s a little function **matsolve(A,b)** that solves for general using the built-in Gaussian Elimination routine **echelon()**, with the addition of a homemade **backsolve()** function. The function in turn relies on a little pivot column detector **pivot()** and my matrix dimension utility **matsize()**. This should include the possibilities of non-square , non-invertible , and treats the case of non-unique solutions in a more or less systematic way.

matsolve(A,b):=block( [AugU], AugU:echelon(addcol(A,b)), backsolve(AugU) ); backsolve(augU):=block( [i,j,m,n,b,x,klist,k,np,nosoln:false], [m,n]:matsize(augU), b:col(augU,n), klist:makelist(concat('%k,i),i,1,n-1), k:0, x:transpose(matrix(klist)), for i:m thru 1 step -1 do ( np:pivot(row(augU,i)), if is(equal(np,n)) then (nosoln:true,return()) else if not(is(equal(np,0))) then (x[np]:b[i], for j:np+1 thru n-1 do x[np]:x[np]-augU[i,j]*x[j]) ), if nosoln then return([]) else return(expand(x)) )$ matsize(A):=[length(A),length(transpose(A))]$ pivot(rr):=block([i,rlen], p:0, rlen:length(transpose(rr)), for i:1 thru rlen do( if is(equal(part(rr,1,i),1)) then (p:i,return())), return(p) )$