Software tools for medium-scale systems of polynomial equations I am attempting to find all real solutions of a system of 12 polynomial equations in 12 unknowns. The equations each have total degree 6 and contain up to 1700 terms. I am only interested in real solutions. The equations were derived as the gradients of a sum-of-squares cost function, which I am attempting to find all global optima of. I believe there are a finite number of real solutions but I have not confirmed this yet. I have floating point coefficients and I'm looking for numerical solutions (as opposed to symbolic solutions).
Which software packages (and which functions specifically) are generally most promising to solve such a problem?
I am aware of various functions in Maple, Matlab, and Mathematica that can solve systems of polynomial equations but there are a large number of options in each software package and I'm interested in advice on where I should be looking first for problems of this scale.
A numerical dump of the cost function is here: 
https://raw.githubusercontent.com/alexflint/polygamy/master/out/epipolar_accel_bezier3/cost.txt
 A: To solve a polynomial system, I would try Bertini which is a homotopy-continuation numerical solver that parallelizes extremely well.  You can also try to attack the optimization problem directly with semi-definite programming as explained by Dima Pasechnik.
A: As far as I can see, you'd like to find a global minimum of $F_0(x)=\sum_{k=1}^{M_0} f_k(x)^2$, where $x=(x_1,\dots,x_{12})$. Equivalently, the problem is to find maximal $\mu$ so that $F(\mu,x):=F_0(x)-\mu\geq 0$ for all $x_1,\dots,x_{12}$.
Now, a sufficient condition for $F(\mu,x)\geq 0$ is that
$F(\mu,x)=\sum_{k=1}^{M} g_k(\mu,x)^2$. The latter condition can be rewritten in the form $F(\mu,x)=\tilde{x}^\top A\tilde{x}\geq 0$ for $\tilde{x}$ being the vector of monomials in $\mu$ and $x$ of degree at most 3, and  $A$ being a symmetric positive semidefinite matrix. The condition $F(\mu,x)=\tilde{x}^\top A\tilde{x}$ is just a bunch of linear equations derived from the coefficients  of the polynomials on the RHS and LHS being equal.
This can be converted into a "semidefinite programming problem" (SDP) to find such a maximal $\mu$, as proposed in Lasserre's  paper; software to do this, and more, is readily available, too; say, YALMIP. The answer is, however, only a lower bound on true value of $\mu$. But often enough it is the exact value what one will get this way.
A: I believe INTLAB, an interval analysis package in MATLAB might be worth considering. It is explained here. Specifically, section 7.4 All Solutions of a Nonlinear System (with an implementation in Appendix A.3) is worth taking a look at.
