You can use the resultant to eliminate a single variable from a system of two polynomials. Unfortunately it seems that these polynomials can't use coefficients from the symbolic ring, so your equations will have to be polynomials over some number field like QQ.

Using your example:

```
sage: var('x, y, z')
sage: eqns = [ x^2 + y^2 == z, x == y]
sage: vars = union(flatten([e.variables() for e in eqns]))
sage: R = PolynomialRing(QQ, vars)
sage: p1, p2 = [R(e.rhs() - e.lhs()) for e in eqns]
sage: SR(p1.resultant(p2, R(z))) == 0
-x + y == 0
sage: SR(p1.resultant(p2, R(x))) == 0
-2*y^2 + z == 0
```

The first two lines are the input as you provided it. The third line collects all the variables occuring in your equations. We use them to construct a multivariate polynomial ring over the rational numbers. We then turn your expressions into polynomials in this ring: each polynomial will be zero iff the corresponding equation holds.

Then we compute the resultant of the two polynomials, with respect to the variable you want to eliminate (`z`

in your question, but `x`

in some other answers). The resultant will be zero exactly if the two polynomials have a common zero. Which in this case means that there is a value for `z`

resp. `x`

which satisfies both equations. So the condition of this polynomial being zero corresponds to the two original equations with that one variable eliminated.

I'm not sure how to generalize this to more than two equations without explicitely solving one equation for one variable.