SymPy Gamma

SymPy:

solve_poly_system ([y**2-x**3+1,y*x],x,y)

$\mathrm{Solve~} \begin{cases} - x^{3} + y^{2} + 1\\x y\end{cases} \mathrm{~for~}\left( x, \ y\right)$

Result:

solve_poly_system ([y**2-x**3+1,y*x],x,y)

See what Wolfram|Alpha has to say.

Experiment with SymPy at SymPy Live.