The problem essentially reduces to finding the roots of a system of polynomials. This article suggests solving the system basically with Newtons method (general

The problem essentially reduces to finding the roots of a system of polynomials. This article suggests solving the system basically with Newtons method (generalized to higher dimensions). But an alternative approach is to use Groebner bases[1] which can give you not just a single (arbitrary) solution (which is what Newtons method produces) but all solutions. SymPy has a good implementation of groebner bases if you want to try it out. [1] http://www.scholarpedia.org/article/Groebner_basis http://www.scholarpedia.org/article/Groebner_basis