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

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!

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}] • Thank you!! I managed to get rid of the holes in the edited solution by increasing the number of datapoints, now it's perfect. Jan 17, 2022 at 17:04
• Still, I'm curious about the original solution. Is there a way to show with Plot3D only the relevant solutions (all variables z1, z2, z3 and z4 being real and non-negative)? Jan 17, 2022 at 17:07