What you can do is choose a monomial ordering on the polynomial ring such that y is the largest variable, so that in the multivariate division algorithm computing the remainder after division by y2−(x3+x) will amount to substituting y2=x3+x:
sage: R.<x,y> = PolynomialRing(GF(43), order='invlex')
sage: (y^2 + x*y + 1).reduce([y^2 - (x^3 + x)])
x*y + x^3 + x + 1
Here you could also choose e.g. R.<y,x> = PolynomialRing(GF(43), order='lex')
.