Polynomial Substitution for only higher terms

edit

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 $y^2 - (x^3 + x)$ will amount to substituting $y^2 = x^3 + 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').

If your polynomial is assumed to be zero, then it is also possible to entirely eliminate $y$ by computing resultant of your polynomial and $y^2-(x^3+x)$ with respect to $y$:

sage: R.<x,y> = PolynomialRing(GF(43))                                                                                                                                                                                 
sage: (y^2+x*y+1).resultant(y^2-(x^3+x),y)                                                                                                                                                                         
x^6 - x^5 + 2*x^4 + x^3 + x^2 + 2*x + 1