rref(a):=block([p,q,k],[p,q]:matrix_size(a),a:echelon(a),
    k:min(p,q),
    for i thru min(p,q) do (if a[i,i]=0 then (k:i-1,return())),
    for i:k thru 2 step -1 do (for j from i-1 thru 1 step -1 do a:rowop(a,j,i,a[j,i])),
    a)$

Degree of polynomial curve
and points to fit

deg:1$
pts:matrix([0,13],[1,13],[2,7],[3,3],[4,5]);

A, b,
check consistency with rref(A|b)

xs:transpose(pts)[1];
ys:transpose(pts)[2];
b:transpose(pts)[2];
np:length(xs);
genmatrix(lambda([i,j],xs[i]^j),np,deg);
A:addcol(col(pts,1)*0+1,%);
addcol(A,ys);
rref(%);

Transpose of A, normal equation At A x = At b,
rref(At A | At b), polynomial curve

transpose(A);
addcol(transpose(A).A,transpose(A).b);
rref(%); 
s:last(transpose(%));
xp:create_list(x^k,k,0,deg);
f(x):=s.xp$
f(x);

Plot the points and the curve

a:first(xs)$
b:last(xs)$
plot2d([f(x),[discrete,xs,ys]],
[x,a,b],[style,lines,points],[point_type,diamond])$