# How do I find all the solutions of three simultaneous equations within a given box?

Sometimes, one needs to find all the solutions of three simultaneous nonlinear equations in three unknowns

\begin{align*}f(x,y,z)&=0\\g(x,y,z)&=0\\h(x,y,z)&=0\end{align*}

within a given cuboidal domain; that is, all triples $(x,y,z)$ satisfying the three equations given above, and within the region defined by $x_{\min}\leq x\leq x_{\max}$, $y_{\min}\leq y\leq y_{\max}$, $z_{\min}\leq z\leq z_{\max}$. (I restrict the discussion here to transcendental equations; algebraic equations are not too problematic for Mathematica (Solve[]/NSolve[], Resultant[], GroebnerBasis[]...))

How can I use Mathematica to find these solutions? FindRoot[] can only find one solution, and you still need an approximate location as a starting point for FindRoot[]. NSolve[] works (sometimes), but it takes long.

-

Stan Wagon's approach of using contour plotting to find solutions to simultaneous equations, as embodied in his function FindAllCrossings2D[], can be generalized to the three-dimensional case. The key is to use ContourPlot3D[] to generate space curves representing the intersection of two of the three functions, as embodied by Maxim Rytin's method (given here). Having done so, one can now use a strategy similar to what was done in FindAllCrossings2D[]; that is, evaluating the third function along the space curve(s), and locating sign changes. The locations of these sign changes can then be used as starting points for FindRoot[].
Here, then, is the routine FindAllCrossings3D[], whose use is completely analogous to its two-dimensional counterpart:
Options[FindAllCrossings3D] =
PerformanceGoal :> $PerformanceGoal, PlotPoints -> Automatic}]]; FindAllCrossings3D[funcs_?VectorQ, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, {z_, zmin_, zmax_}, opts___] := Module[{contourData, seeds, tt, fz = Compile[{x, y, z}, Evaluate[funcs[[3]]]]}, contourData = Cases[Normal[ContourPlot3D[ Evaluate[Most[funcs]], {x, xmin, xmax}, {y, ymin, ymax}, {z, zmin, zmax}, BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {}}, ContourStyle -> None, Mesh -> None, Method -> Automatic, Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings3D]], Options[ContourPlot3D]]]]], Line[l_] :> l, Infinity]; seeds = Flatten[Pick[Rest[#], Most[#] Rest[#] &@Sign[Apply[fz, #, 2]], -1] & /@ contourData, 1]; If[seeds === {}, seeds, Select[Union[Map[{x, y, z} /. FindRoot[funcs, Transpose[{{x, y, z}, #}], Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings3D]], Options[FindRoot]]]] &, seeds]], (xmin < #[[1]] < xmax && ymin < #[[2]] < ymax && zmin < #[[3]] < zmax) &]]]  As an example of how to use FindAllCrossings3D[]: sols = FindAllCrossings3D[ {Sin[x + y] Sin[y - z], Cos[x] Cos[y] - Sin[z], x^2 + y^2 + z^2 - 9}, {x, -4, 4}, {y, -4, 4}, {z, -4, 4}] {{-2.80293, -0.756176, -0.756176}, {-2.78082, -0.360773, -1.06625}, {-2.11276, 2.11276, 0.269309}, {-1.14056, -0.395145, 2.74645}, {-1.14056, 2.74645, -0.395145}, {-0.883563, 0.883563, 2.72739}, {-0.360773, -2.78082, -1.06625}, {0.360773, 2.78082, -1.06625}, {0.883563, -0.883563, 2.72739}, {1.14056, -0.395145, 2.74645}, {1.14056, 2.74645, -0.395145}, {2.11276, -2.11276, 0.269309}, {2.78082, 0.360773, -1.06625}, {2.80293, -0.756176, -0.756176}}  The routine found$14$solutions. To visualize the solutions, we can do the following: l1 = Cases[Normal[ContourPlot3D[{Sin[x + y] Sin[y - z], Cos[x] Cos[y] - Sin[z]}, {x, -4, 4}, {y, -4, 4}, {z, -4, 4}, BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {}}, ContourStyle -> None, Mesh -> None]], Line[l_] :> l, Infinity]; Graphics3D[{Line[l1], Sphere[{0, 0, 0}, 3], Sphere[sols, 1/10]}, Axes -> Automatic]  where we used small spheres to mark the intersections of the space curves formed by the intersection of$\sin(x+y)\sin(y-z)=0$and$\cos\,x\cos\,y=\sin\,z$, and the sphere$x^2+y^2+z^2=9\$.