gap> R:=PolynomialRing(Rationals,["x","y"]);
Rationals[x,y]
gap> AssignGeneratorVariables(R);
#I  Assigned the global variables [ x, y ]
gap> gens:=[x*y^3-x^2,x^3*y^2-y];
[ x*y^3-x^2, x^3*y^2-y ]
gap> order:=MonomialGrlexOrdering();
MonomialGrlexOrdering()
gap> bas:=ReducedGroebnerBasis(gens,order);
[ y^4-x*y, x*y^3-x^2, x^4-y^2, x^3*y^2-y ]

# the possible remainders: x,y- degree is less than 4 and cannot reduce:
gap> remainders:=List(Cartesian([0..3],[0..3]),a->x^a[1]*y^a[2]);
[ 1, y, y^2, y^3, x, x*y, x*y^2, x*y^3, x^2, x^2*y, x^2*y^2, x^2*y^3, x^3, 
  x^3*y, x^3*y^2, x^3*y^3 ]
                                                                            ^
gap> remainders:=Filtered(remainders,
       a->a=PolynomialReducedRemainder(a,bas,order));
[ 1, y, y^2, y^3, x, x*y, x*y^2, x^2, x^2*y, x^2*y^2, x^3, x^3*y ]

# calculate (reduced) images when multiplying with x,y

gap> ximgs:=List(remainders,
> a->PolynomialReducedRemainder(a*x,bas,order));
[ x, x*y, x*y^2, x^2, x^2, x^2*y, x^2*y^2, x^3, x^3*y, y, y^2, y^3 ]
gap> yimgs:=List(remainders,
> a->PolynomialReducedRemainder(a*y,bas,order));
[ y, y^2, y^3, x*y, x*y, x*y^2, x^2, x^2*y, x^2*y^2, x^3, x^3*y, y ]

# to build matrices we need to write in terms of basis:

gap> TermsOfPolynomial:=function(pol)
> local t,c,l;
> l:=[];
> while not IsZero(pol) do
> t:=LeadingMonomialOfPolynomial(pol,order);
> c:=LeadingCoefficientOfPolynomial(pol,order);
> pol:=pol-c*t;
> Add(l,[c,t]);
> od;
> return l;
> end;
function( pol ) ... end
gap> TermsOfPolynomial(gens[1]);
[ [ 1, x*y^3 ], [ -1, x^2 ] ]

# and build a vector

gap> decompose:=function(pol)
> local c,i;
> c:=ListWithIdenticalEntries(Length(remainders),0);
> for i in TermsOfPolynomial(pol) do
> c[Position(remainders,i[2])]:=i[1];
> od;
> return c;
> end;
function( pol ) ... end

# now write images as rows in terms of basis
gap> xmat:=List(ximgs,decompose);
[ [ 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ], 
  [ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 ] ]

gap> ymat:=List(yimgs,decompose);
[ [ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 ], 
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ], 
  [ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] ]