# Plotting the number of fixed points of a system of two nonlinear differential equations

I am working with a complex set of two differential equations with many parameters. A simplified minimal working example would be the following system:

$$x'[t]=x^2[t]+y[t]+b$$ $$y'[t]=x[t]+y^2[t]-a$$

,where $a\in [0,1]$ and $b\in [0,1]$ are two parameters.

Now, depending on parameters the system may have 0,1 or 2 fixed points. What I want to do is Regionplot/Contourplot the parameterspace $[0,1]$x$[0,1]$ and shade regions/plot the boundary of regions for which there exist a different number of equilibria.

What I started with is plotting the nullclines as well as computing the fixed points of the system. I could just count the number of fixed points but I dont know how much that will help.

a = 0.5;
b = 0.5;
eq1=x^2 + y + b;
eq2=x + y^2 - a;
ContourPlot[{eq1 == 0, eq2 == 0}, {x, -3, 3}, {y, -3,
3}, ContourStyle -> {Red, Blue}, ImageSize -> Scaled[0.3],
PlotLegends ->
LineLegend[{"\!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]\)x=0",
"\!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]\)y=0"}],
FrameLabel -> {x, y}]

a = NSolve[{eq1 == 0 && eq2 == 0 }, {x, y}, Reals]
Length[a]


Ideally, I would like to end up with a plot similar to the following, where either regions are shaded or the boundary-line is plotted: Edit: I want to emphasize that I am looking for numerical methods for a more complex problem.

Here is my current code, which works but still has 2 major issues:

1.I was only able to ListPlot3D the list and I would like to have a 2D projection of the data.

2.Empty sets, i.e. parameter configurations for which there exist no fixed point are still counted. I would like to get rid of those entries, while still preserving the value 0 in the plot.

eq1 = x^2 + y + b;
eq2 = x + y^2 - a;

c = Table[{a, b, x} /.
NSolve[{eq1 == 0 && eq2 == 0 }, {x, y}, Reals], {a, 0, 1,
0.005}, {b, 0, 1, 0.005}];

d = Flatten[Map[Length, c, {2}]];
e = Table[{a, b}, {a, 0, 1, 0.005}, {b, 0, 1, 0.005}];
f = Flatten[e, 1];
g = Flatten[Riffle[f, d]];
par = Partition[g, 3];

ListPointPlot3D[par, ColorFunction -> "Rainbow", Boxed -> False,
Axes -> {True, True, False}, ViewPoint -> {0, 0, 250},
PlotLegends -> SwatchLegend[{Red, Blue}, {0, 2}]] • Your a has two meanings, which is usually a bad idea. – corey979 Jan 24 '17 at 16:10

Update

Inspired by Chris K and showing larger parameter range (full range of real critical points: 0,1,2,3,4:

fun[a_, b_] :=
Length[{x, y} /.
NSolve[{f[a, b, x, y] == {0, 0}}, {x, y}, Reals] /. {x, y} :> {}];
g[a_, b_] := (27 + 256 a^3 + 288 a b - 256 a^2 b^2 - 256 b^3);
h[a_, b_] := If[g[a, b] < 0, fun[a, b], 2];
lcp = ContourPlot[h[a, b], {a, -8, 8}, {b, -8, 8},
PlotLegends -> SwatchLegend[{Red, Green, Blue}, {0, 2, 4}]];
DynamicModule[{p = {1/2, 1/2}},
Row[{LocatorPane[Dynamic[p],
Show[lcp, Graphics[{PointSize[0.03], Yellow, Dynamic[Point[p]]}],
ImageSize -> 300]],
Dynamic@ContourPlot[
Evaluate@Thread[(f[##, x, y] & @@ p) == 0], {x, -4, 4}, {y, -4,
4}, ContourStyle -> {Red, Blue}, ImageSize -> 300]}]] For these quadratic relations it is relatively straight forward to use the discriminant, e.g.:

f[a_, b_, x_, y_] := {x^2 + y + b, x + y^2 - a}
DynamicModule[{p = {1/2, 1/2}},
Row[{LocatorPane[Dynamic[p],
ContourPlot[
27 + 256 a^3 + 288 a b - 256 a^2 b^2 - 256 b^3, {a, 0, 1}, {b, 0,
1}, Contours -> {{0}}, ContourShading -> {Red, Blue},
ContourStyle -> Directive[Yellow, Thickness[0.02]],
Epilog -> {PointSize[0.03], Yellow, Dynamic[Point[p]]},
ImageSize -> 300]],
Dynamic@
ContourPlot[
Evaluate@Thread[(f[##, x, y] & @@ p) == 0], {x, -2, 2}, {y, -2,
2}, ContourStyle -> {Red, Blue}, ImageSize -> 300]
}]
]


An example of no real zeroes: "Approximately" 1: Two real zeros: • Thank you very much. This solves the problem for the minimal working example in the description. However, the problem I am dealing with contains polynomials of 5th order and hence I was not able to derive an analytical solution. It was kind of my bad in posting this overly simplified working example above. Could anyone try to get this working numerically? – ParetoWilli Jan 25 '17 at 8:50
• @ParetoWilli yes, I'd also like to see a purely numerical approach! – Chris K Jan 25 '17 at 13:23

For this particular example you can find the bifurcations analytically by also setting the slopes of the isoclines to be equal.

eq1 := x^2 + y + b;
eq2 := x + y^2 - a;
bif = Solve[{eq1 == 0, eq2 == 0, D[eq1, x]/D[eq1, y] == D[eq2, x]/D[eq2, y]}, {x, y, a, b}]

(* {{y -> 1/(4 x), a -> (1 + 16 x^3)/(16 x^2), b -> (-1 - 4 x^3)/(4 x)}} *)

ParametricPlot[{a, b}/. bif, {x, -10, 10}, PlotRange -> {{-4, 8}, {-8, 4}}, Frame -> True, Axes -> False] I couldn't resist looking at a larger range of a and b, because there is another bifurcation between 2 and 4 equilibria. I don't have a way to automatically label the number of equilibria, but in the upper left there are 0, in the middle there are 2, and in the lower right there are 4.

Taking a look at the phase planes at three bifurcation points where a=2:

bif2 = NSolve[{eq1 == 0, eq2 == 0, D[eq1, x]/D[eq1, y] == D[eq2, x]/D[eq2, y], a == 2}, {x, y, a, b}]

(* {{a -> 2., x -> 1.98412, y -> 0.126, b -> -4.06275}, {a -> 2., x -> -0.169722, y -> -1.473, b -> 1.44419}, {a -> 2., x -> 0.185598, y -> 1.347, b -> -1.38144}} *)

GraphicsRow[Table[{a, b} = {a, b} /. bif2[[i]]; ContourPlot[{eq1 == 0, eq2 == 0}, {x, -3, 3}, {y, -3, 3},
ContourStyle -> {Red, Blue}], {i, 3}]] • best answer. +1 an ashamed of mine though if I extend my parameter range my answer agrees with yours though my 'discriminant' does not allow different shading for 3 and 4 root case... – ubpdqn Jan 25 '17 at 4:16
• @ubpdqn There is still plenty of room for improvement to really answer the question: 1) counting # of equilibria, 2) shading regions, 3) numerical methods for more complex models. Hopefully someone takes a stab at it. – Chris K Jan 25 '17 at 4:20