Sage (using the interface to Singular) can solve multivariate polynomial equations in some situations (they assume that the solutions form a zero-dimensional variety) using Gröbner bases. Here is a simple example:
sage: R = PolynomialRing(QQ, 2, 'ab', order='lp') sage: a,b = R.gens() sage: I = (a^2-b^2-3, a-2*b)*R sage: B = I.groebner_basis(); B [b^2 - 1, a - 2*b]
See About this document... for information on suggesting changes.