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

Question

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}]]

Answer 1

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}, 
   ContourShading -> {Red, Green, Blue}, 
   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]}]]

Original Answer

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:

Answer 2

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}]]