Return to Question

Edit 4: Just for fun:

addConstraints[foc_] := {Sequence @@ foc, x >= 1, y >= 1, z >= 1}
Reduce[addConstraints@Thread[foc == 0], {x, y, z}, Backsubstitution -> True] // N // Timing
{6.87911, x == 6.18083 && y == 4.58414 && z == 3.14352}

gb := GroebnerBasis[foc, {x, y, z}];
NSolve[addConstraints[FactorList[#][[-1, 1]] & /@ gb], Reals] // Timing
{2.97643, {{x -> 6.18083, y -> 4.58414, z -> 3.14352}}}

Solve[addConstraints@Thread[foc == 0]] // N // Timing
{6.84882, {{x -> 6.18083, y -> 4.58414, z -> 3.14352}}}

Edit 3: belisarius has pointed out that foc/.{x->0,y->0} = {0,0,0}. It is also the case that foc/.{x->0,z->0} = {0,0,0} and foc/.{y->0,z->0} = {0,0,0}. I cannot trivially delete these roots from the original equations, because the original equations do not factor so. However, these roots do seem to present the problem when using NSolve after GroebnerBasis as is evident from the error message. So, I've tried dropping all the small common factors from the Groebner basis.

gb = GroebnerBasis[foc,{x,y,z}];
NSolve[FactorList[#][[-1,1]]&/@gb,Reals]

gives the set of the remaining solutions, with the one I was looking for being among them.

So, it seems that NSolve gets caught in the positive dimensionality part of the solution, even if explicitly told the region (x>1,y>1,z>1), where that positive dimensionality is not present. When run after GroebnerBasis, NSolve at least admits it, but when run on the original equations, it seems to hide the fact.

(I accept belisarius' answer as the comment below helped me to at least have an idea of what might be going wrong.)

Edit 2: As is noted in the comments, the system above is huge, and it takes too long to try different solution methods. Interestingly, while working on the initial problem, I came across a much simpler system, but with the same NSolve issue. Here are the equations:

foc = {-x^4 y^2+2 x^3 y^3+x^2 y^4-2 x^4 y z+4 x^3 y^2 z-x^4 y^2 z+10 x^2 y^3 z-2 x^3 y^3 z+4 x y^4 z-x^2 y^4 z-2 x^4 z^2+8 x^2 y^2 z^2+8 x y^3 z^2+2 y^4 z^2,x^4 y^4+2 x^3 y^5-x^2 y^6+8 x^4 y^3 z+18 x^3 y^4 z-2 x^4 y^4 z+4 x^2 y^5 z-4 x^3 y^5 z-2 x y^6 z-2 x^2 y^6 z+16 x^4 y^2 z^2+44 x^3 y^3 z^2-4 x^4 y^3 z^2+28 x^2 y^4 z^2-10 x^3 y^4 z^2-8 x^2 y^5 z^2-2 y^6 z^2-2 x y^6 z^2+12 x^4 y z^3+40 x^3 y^2 z^3-2 x^4 y^2 z^3+38 x^2 y^3 z^3-6 x^3 y^3 z^3+8 x y^4 z^3-7 x^2 y^4 z^3-2 y^5 z^3-4 x y^5 z^3-y^6 z^3+3 x^4 z^4+12 x^3 y z^4+14 x^2 y^2 z^4+4 x y^3 z^4-y^4 z^4,y (2 x^2 y^3+8 x^2 y^2 z+4 x y^3 z+10 x^2 y z^2+10 x y^2 z^2-2 x^2 y^2 z^2+y^3 z^2-x y^3 z^2+4 x^2 z^3+6 x y z^3-4 x^2 y z^3+2 y^2 z^3-4 x y^2 z^3-2 x^2 z^4-3 x y z^4-y^2 z^4)}

Please note that these are not symmetric anymore. NSolve[foc,Reals] gives {} without any errors. So does NSolve[Join[foc, {x >= 1, y >= 1, z >= 1}], Reals] (the constraints give the domain I am interested in). Running GroebnerBasis before NSolve still produces no results but the following error message appears from NSolve:

Infinite solution set has dimension at least 1.
Returning intersection of solutions with
(171802 x)/178835-(113492 y)/178835-(121484 z)/178835 == 1.

The best part of it, the solution exists as can be easily seen from the following contour plot (and by playing with the upper limits it seems the solution is unique):

ContourPlot3D[Evaluate@Thread@(foc == 0), {x, 1, 10}, {y, 1, 10}, {z, 1, 10}]

Solve, given the constraints on the domain, can find that unique solution. Somehow, NSolve fails in this case...