1. Harvesting the Seeds
Perhaps the first thing one wants to know about an autonomous system dx/dt =
f(x, y), dy/dt = g(x, y) is where the equilibrium points are. For such a point,
dx/dt and dy/dt are both zero, so the problem is the same as solving a pair of
simultaneous equations. Standard root-finding methods can do this almost
instantaneously provided they are given a starting value not too far away from
the root. Here is an example. First we define the system that will be our main
example, and then we appeal to FindRoot.
In[]:=

The root is returned in the form of a substitution rule. We can define rx and
ry to be the x and y coordinates of the root, and then we plug into the given
functions as a check (% refers to the preceding output).
In[]:=

f[x_, y_] := x - y^2 Cos[y]
g[x_, y_] := -y + x Sin[x]
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, -6}, {y, 3}]
{rx, ry} = {x, y} /. %
{f[rx, ry], g[rx, ry]}
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

So far, so good. But our problem here is a little different in that we want to
know all the equilibrium points in a given rectangle. Here is where the power of
Mathematica , and its openness in allowing us access to the data in a
contour plot comes into play. Consider the following algorithm:
1. Generate the curve(s) representing f = 0 by using a contour plot for the 0
contour.
2. Evaluate g along the curves, and check for sign changes of g.
3. Since a sign change of g along f = 0 means that we are quite close to a
root,
we can use the gathered values as seeds to FindRoot.
In fact, this procedure works very well and is not too difficult to program.
The function Equilibria which accompanies this article implements it in detail.
Since the main outline is so simple, we illustrate the programming methodology
by following through an example. Graphics[ContourPlot[***]] returns a
complicated graphics object that has imbedded in it the several Line objects
making up the contours. Cases[***, _Line, Infinity] picks them out no matter
where they are. First then strips off the Line header because we are really
interested in these objects as lists. We reintroudce the Line headers in order
to confirm visually that this scheme has worked.
In[]:=

f[x_, y_] := x - y^2 Cos[y]
contourData = Map[First,
Cases[Graphics[
ContourPlot[f[x, y], {x, -10, 10}, {y, -10, 10},
Contours -> {0}, PlotPoints -> 50]],
_Line, Infinity]];
Show[Graphics[Map[Line, contourData]]];
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

Good. We now have the five curves as lists in contourData. Let us now focus on
the central curve, given by contourData[[5]]. We apply g and look at the signs
of the result.
In[]:=

f[x_, y_] := x - y^2 Cos[y]
g[x_, y_] := -y + x Sin[x]
contourData = Map[First,
Cases[Graphics[
ContourPlot[f[x, y], {x, -10, 10}, {y, -10, 10},
Contours -> {0}, PlotPoints -> 50]],
_Line, Infinity]];
signs = Sign[Apply[g, contourData[[5]], 2]] // Short
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

Multiplying this vector by a shift of it causes the sign-changes to show up as
-1s. There are six of them.
In[]:=

f[x_, y_] := x - y^2 Cos[y]
g[x_, y_] := -y + x Sin[x]
contourData = Map[First,
Cases[Graphics[
ContourPlot[f[x, y], {x, -10, 10}, {y, -10, 10},
Contours -> {0}, PlotPoints -> 50]],
_Line, Infinity]];
signs = Sign[Apply[g, contourData[[5]], 2]] // Short
Rest[%] * Rest[RotateRight[%]] // Short
Position[%, -1]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

We now find the places in the plane where these sign-changes occurred.
In[]:=

f[x_, y_] := x - y^2 Cos[y]
g[x_, y_] := -y + x Sin[x]
contourData = Map[First,
Cases[Graphics[
ContourPlot[f[x, y], {x, -10, 10}, {y, -10, 10},
Contours -> {0}, PlotPoints -> 50]],
_Line, Infinity]];
signs = Sign[Apply[g, contourData[[5]], 2]] // Short
Rest[%] * Rest[RotateRight[%]] // Short
Position[%, -1]
seeds = contourData[[5, Flatten[%]]]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

This gives us six approximations. We define a function, root, that feeds a
single seed to the rootfinder; mapping this function to the set of seeds
produces six fairly exact roots.
In[]:=

f[x_, y_] := x - y^2 Cos[y]
g[x_, y_] := -y + x Sin[x]
contourData = Map[First,
Cases[Graphics[
ContourPlot[f[x, y], {x, -10, 10}, {y, -10, 10},
Contours -> {0}, PlotPoints -> 50]],
_Line, Infinity]];
signs = Sign[Apply[g, contourData[[5]], 2]] // Short
Rest[%] * Rest[RotateRight[%]] // Short
Position[%, -1]
seeds = contourData[[5, Flatten[%]]]
root[{seedx_, seedy_}] :=
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, seedx}, {y, seedy}]
Map[root, seeds]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

The following figure, which shows the contours together with circles around the
data we just found, shows that the approach has worked beautifully and we have
found all six solutions on the central contour.
In[]:=

f[x_, y_] := x - y^2 Cos[y]
g[x_, y_] := -y + x Sin[x]
contourData = Map[First,
Cases[Graphics[
ContourPlot[f[x, y], {x, -10, 10}, {y, -10, 10},
Contours -> {0}, PlotPoints -> 50]],
_Line, Infinity]];
signs = Sign[Apply[g, contourData[[5]], 2]] // Short
Rest[%] * Rest[RotateRight[%]] // Short
Position[%, -1]
seeds = contourData[[5, Flatten[%]]]
root[{seedx_, seedy_}] :=
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, seedx}, {y, seedy}]
Map[root, seeds]
Show[ContourPlot[#[x, y], {x, -10, 10}, {y, -10, 10},
Contours->{0}, ContourShading->False,
PlotPoints ->80] & /@ {f, g},
Graphics[{Circle[{x, y}, 0.32] /. %} ]];
Show[ContourPlot[#[x, y], {x, -10, 10}, {y, -10, 10},
Contours->{0}, ContourShading->False,
DisplayFunction->Identity, PlotPoints->100,
ContourStyle->AbsoluteThickness[.6]] & /@ {f, g},
Graphics[{AbsoluteThickness[.75],
Circle[{x, y}, 0.32] /. %%} ],
DisplayFunction->$DisplayFunction,
PlotRange->All, AspectRatio->Automatic];
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

As mentioned, Equilibria implements this in general. It finds all 16 equilibria
points for our example in short order.
In[]:=

<<Local`NullclinePlot`
f[x_, y_] := x - y^2 Cos[y]
g[x_, y_] := -y + x Sin[x]
Equilibria[{f[x, y], g[x, y]}, {x, -10, 10}, {y, -10, 10}]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above

There are some options to Equilibria. Setting SeedsOnly to True causes just the
seeds to be returned. The EquilibriaPlotPoints option controls the number of
plot points in the contour plots that get the algorithm started; for difficult
examples this might have to be increased from the default of 49. And there is a
ShowEigenvalues option that produces a table of the equilibrium points together
with the eigenvalues of the Jacobian at the points.
EXERCISE. Modify these ideas to get a routine that finds all complex roots of a
complex function f(z) in a given rectangle. Just apply the same method, looking
at f's imaginary part along a contour plot of the real part. FindRoot, does in
fact work on complex functions, provided I appears in either the function or the
initial value.
Up to New Visualization Ideas for Differential Equations