Show solutions of non-linear equations with multiple variables in 3D space

Question

Consider a (non-linear) system of equations with two parameters x and y and multiple variables z1, z2, z3, z4

 eqs = {
   -0.001 z1 + (0.001 + 0.1 z1) (1 - z1 - z2 - z3 - z4) - 0.06 z1 (z3 + z4) == 0,
   0.001 z1 - 0.011 z2 + (0.05 + x (1 - y)) z3 - 0.6 z2 (z3 + z4) == 0,
   (1 - z1 - z2 - z3 - z4) (0.001 + 0.6 z3 + 0.6 z4) - (0.05 + x) z4 == 0,
   0.001 z2 - (0.06 + x (1 - y)) z3 + 0.06 z1 (z3 + z4) + 0.6 z2 (z3 + z4) == 0
   };

The system has as many equations as variables (4) and can be solved given values for the parameters

NSolve[Join[eqs /. {x -> .6, y -> .5}, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1, z2, z3, z4}]

Out: {{z1 -> 0.897951, z2 -> 0.0915341, z3 -> 0.000372715, z4 -> 0.0000192353}}

Sometimes it has more than one solution

NSolve[Join[eqs /. {x -> .7, y -> .8}, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1, z2, z3, z4}]

Out:
{
 {z1 -> 0.894114, z2 -> 0.0942418, z3 -> 0.00107801, z4 -> 0.0000233464},
 {z1 -> 0.274034, z2 -> 0.256676, z3 -> 0.249387, z4 -> 0.0417326},
 {z1 -> 0.636264, z2 -> 0.250085, z3 -> 0.0626664, z4 -> 0.00259097}
}

Now I want to show the solutions for z1 in the x,y plane. The following works ok

list = Flatten[Table[{x, y, z1 /. #} & /@ NSolve[Join[eqs, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1, z2, z3, z4}], {x, 0, 1, .025}, {y, 0, 1, .025}], 2];
ListPointPlot3D[list, ColorFunction -> "Rainbow", AxesLabel -> {"x", "y", "z1"}]

But I find it a bit difficult to interpret. I was hoping to show the solution as continuous planes, as in the help file for ContourPlot3D

However I don't know how to use ContourPlot3D because of the other variables in the system z2, z3, z4. My first idea was to reduce the system to a single equation with just z1. That didn't work though because the system is non-linear and Solve didn't find a solution.

Any ideas?

Thanks!

Answer 1

Edit

The surface still have some holes.

eqs = {-0.001 z1 + (0.001 + 0.1 z1) (1 - z1 - z2 - z3 - z4) - 
     0.06 z1 (z3 + z4) == 0, 
   0.001 z1 - 0.011 z2 + (0.05 + x (1 - y)) z3 - 0.6 z2 (z3 + z4) == 
    0, (1 - z1 - z2 - z3 - z4) (0.001 + 0.6 z3 + 0.6 z4) - (0.05 + 
        x) z4 == 0, 
   0.001 z2 - (0.06 + x (1 - y)) z3 + 0.06 z1 (z3 + z4) + 
     0.6 z2 (z3 + z4) == 0};
eqs = Rationalize[eqs, 0];
data = Flatten[
   Table[{x, y, z1} /. 
     NSolve[Join[
       eqs, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1}, {z2, z3, 
       z4}], {x, 0, 1, .01}, {y, 0, 1, .01}], 2];
surf = ListSurfacePlot3D[data, MaxPlotPoints -> 60, 
  ColorFunction -> "Rainbow", PerformanceGoal -> "Quality", 
  ViewPoint -> {-2.60, -1.74, 1.27}]

Original

The solution have several branch.

eqs = {-0.001 z1 + (0.001 + 0.1 z1) (1 - z1 - z2 - z3 - z4) - 
     0.06 z1 (z3 + z4) == 0, 
   0.001 z1 - 0.011 z2 + (0.05 + x (1 - y)) z3 - 0.6 z2 (z3 + z4) == 
    0, (1 - z1 - z2 - z3 - z4) (0.001 + 0.6 z3 + 0.6 z4) - (0.05 + 
        x) z4 == 0, 
   0.001 z2 - (0.06 + x (1 - y)) z3 + 0.06 z1 (z3 + z4) + 
     0.6 z2 (z3 + z4) == 0};
eqs = Rationalize[eqs, 0];
sol = Solve[eqs, {z1}, {z2, z3, z4}];
Plot3D[z1 /. sol// Evaluate, {x, 0, 1}, {y, 0, 1},MaxRecursion -> 4, ColorFunction -> "Rainbow", 
 AxesLabel -> {"x", "y", "z1"}, BoxRatios -> {1, 1, 1}]