My problem is to find the minimum of a function $f({\bf x})$ given a set of constraints $c({\bf x})={\bf 0}$. I could get the equations using lagrange formulation: $$ \triangledown_{({\bf x,\lambda})}(f+ {\bf \lambda} c({\bf x})) = {\bf 0}$$ Depending on the parametrization of my function $f$, the equations can be non-algebraic (invloving angles and their trigonometric functions simultaneously) or algebraic (using unit quaternions at a compromise of increasing the number of equations and variables). Without going into much details of the physical problem, I am solving it by considering a number of finite elements (say $N$). The higher the value of $N$, the better the modeling of the problem, the larger the computational time since the number of equations and variables increase linearly with $N$.
This is a standard optimization problem which can be solved using Matlab, Maple, etc. in seconds for up to $N=100$ given an intial guess of the solution. But, I get only one solution and it could be a local minima, thus, the assumed initial value plays an important role. I don't want to use any intial guess. So, I tried using Bertini (which uses homotopy continuation) to obtain all the solutions for $N=3$, where my number of equations and variables are already 18. It takes a while (about an hour) and gives me all the solutions. But it fails for $N>3$. I have looked into interval analysis methods but again I have to specify intervals for my huge number of variables which I wouldn't know.
I am looking for any suggestions on other algorithms or softwares that can numerically handle a huge number of small degree multi-variate equations with a huge number of variables in a considerable amount of time if at all feasible. The degrees of my equations are 3 at most and I don't think there are any positive dimensional components in the solutions.