ring r=0,(x(1..4)),dp;    //number of variables must be the same as the number of columns of the matrix for which we want the reduced row echelon form

matrix m[4][4]=16,2,3,13,5,11,10,8,9,7,6,12,4,14,15,1; //enter an arbitrary matrix here

print("m=");print(m);

int nr=nrows(m);
int nc=ncols(m);
int j;
int k;
for(j=1;j<=nr;j=j+1){for(k=1;k<=nc;k=k+1){m[j,k]=m[j,k]*x(k);}};
kill j;
kill k;

matrix ONES[nc][1];
for(int j=1;j<=nc;j=j+1){ONES[j,1]=1;};
kill j;

ideal i=m*ONES;
option(redSB);
ideal g=std(i);
int L=size(g);
for(int j=1;j<=L;j=j+1) {g[j]=g[j]/leadcoef(g[j]);} // set the coefficients of the leading terms to unity
kill j;

matrix RREF[L][nc];
int k;
int j;
for (j=L;j>0;j=j-1){for (k=1;k<=nc;k=k+1){RREF[j,k]=jet(g[L+1-j]/x(k),0);}};

print("RREF=");print(RREF);print("RREF is the Reduced Row Echelon Form of the matrix m")

ring R=0,x,dp;
matrix m[4][3]=1,2,3,4,5,6,7,8,9,10,11,12;
list L=ludecomp (m);
print(L[3]);
1,2, 3,
0,-3,-6,
0,0, 0,
0,0, 0