Plotting solutions of ODE (with arrows!)

By Gianluca Gorni, Department of Mathematics and Computer Science, University of Udine, Italy. Latest revision: June 7, 2003.

Description: The simplest functions exported by this package, OrientedParametricPlot and NDSolvePlot, place arrow heads along plane parametric curves, or along solutions to differential equations, to show how they are oriented. The positioning is made with respect to arc length.
For phase portraits of dynamical systems, where many curves are to be displayed together, an  algorithm is set up to automatically place arrows taking into account only the parts of the curves that are within the viewing area, with the number of arrows on each curve increasing with the length. This automatic arrow-placing is used in the functions HamiltonianPlot, PhasePlot, LinearPhasePlot.
The functions PoincareMap and PoincareMapPlot are defined to help visualize the Poincaré map of a two-dimensional system.To assist in finding fixed points of the Poincaré map the following functions are defined: PoincareMapFixedPointDensityPlot, PoincareMapFixedPointPlot3D and FindPoincareMapFixedPoint.

I am aware of two other packages freely available over the internet whose functionality partially overlaps with the present one. They are DrawingArrows by David Park and dynapack by Phillip Kent and Phil Ramsden.

Instructions: load the package by clicking on the button below, then evaluate the examples. This notebook is not meant for absolute beginners of Mathematica. It assumes some familiarity with notebooks and with the syntax of basic standard Mathematica commands such as Plot and ParametricPlot.

Button[
      StyleForm["Click here to load package!", 
        FontColor -> RGBColor[1, 0, 0], FontSize -> 14], 
      ButtonFunction :> 
        FrontEndExecute[FrontEndToken["EvaluateInitialization"]], 
      ButtonEvaluator -> None, Active -> True]

"Click here to load package!"

Mathematica package

This is the program. It is meant to be read only by Mathematica programmers.

Preliminaries and usage documentation

Off[General::"spell"]

Off[General::"spell1"]

\(\(BeginPackage["ODEPlot`", "Utilities`FilterOptions`", "Graphics`Arrow`", "Graphics`PlotField`"];\)\)

Unprotect["`*"];
ClearAll["`*"];

\(\(LinearPhasePlot::usage = 
"LinearPhasePlot[matrix,options] plots a phase portrait of the system of two differential equations {x'[t],y'[t]}==mat.{x[t],y[t]}. The matrix must be numeric and exact. The curves are drawn from the exact symbolic solutions of the system. Options are NotableLineStyle, NumberOfTrajectories, WriteReport, ShowField, CurvesWithArrows, ArrowDensity. LinearPhasePlot also accepts all options for Graphics`Arrow` and Graphics. Example: LinearPhasePlot[{{1,2}},{-2,1}}].";\)\)

NotableLineStyle::usage = "NotableLineStyle is an option for \
LinearPhasePlot specifying how to highlight the straight lines that are of \
special significance for the system (either eigenspaces or principal \
directions). The default value is \
NotableLineStyle->RGBColor[1,0,0].";

NumberOfTrajectories::usage = "NumberOfTrajectories is an option \
for LinearPhasePlot that helps control how many trajectories to (attempt to) \
draw, approximately. The default value is NotableLineStyle->5.";

WriteReport::usage = "WriteReport is an option for LinearPhasePlot \
specifying whether to write the classification of the given system. The \
default value is WriteReport->True.";

ShowField::usage = "ShowField is an option for PhasePlot, \
LinearPhasePlot and HamiltonianPlot specifying whether to draw the direction \
field of the system (by means of the package Graphics`PlotField`) together \
with the trajectories. The default value is ShowField->False.";

CurvesWithArrows::usage = "CurvesWithArrows is an option for \
LinearPhasePlot and PhasePlot, specifying whether to draw arrows on the \
trajectories. The default value is CurvesWithArrows->True. With \
CurvesWithArrows->False the plot takes less to compute.";

ArrowDensity::usage = "ArrowDensity is an option for \
LinearPhasePlot, HamiltonianPlot and and PhasePlot, influencing how many \
arrows to draw along the lines. The default value is ArrowDensity->1. \
ArrowDensity should be set to a positive number. Higher values should yield \
more arrows. You cannot say for sure beforehands how many, because arrow \
placement takes into account the length of the part of the curves that lies \
inside the display area, ignoring the rest.";

\(\(OrientedParametricPlot::usage = "OrientedParametricPlot[{x[t],y[t]},{t,tMin,tMax},opts] makes a ParametricPlot and then adds arrows along the curve. Specific options are HowManyArrows, ArrowPositions, ReverseDirectionQ. It also accepts options from ParametricPlot and Graphics`Arrow`. HeadScaling is set to Absolute by default.";\)\)

HowManyArrows::usage = "HowManyArrows is an option for \
OrientedParametricPlot that directs how many arrows to add along the curve. \
The default is HowManyArrows->1.";

ArrowPositions::usage = "ArrowPositions->{p1,p2,...} is an option \
for OrientedParametricPlot that directs where to draw the arrows heads along \
the curve. The values p1,p2,... are supposed to be numbers between 0 and 1 \
(otherwise they are normalized automatically), and refer to the curve \
parameterized by arc length. The default is ArrowPositions->Automatic, under \
which the lengths of the arcs between successive arrow heads are the same, \
and double the arcs at either end. This arrangement is best for closed cuves.\
";

ReverseDirectionQ::usage = "ReverseDirectionQ is an option for \
OrientedParametricPlot that directs whether to draw the arrows in the \
direction of the given parameterization, or the opposite one. The default is \
ReverseDirectionQ->False.";

TrigInterpolation::usage = "TrigInterpolation[data,x] is ";

EvaluateOnCoefficients::usage = "EvaluateOnCoefficients is \
";

Period::usage = "Period is ";

\(\(TrigInterpolationPlot::usage = "TrigInterpolationPlot[data,xmin,xmax] is ";\)\)

ShowListPlot::usage = "ShowListPlot is ";

ListPlotStyle::usage = "ListPlotStyle is ";

ShowFormula::usage = "ShowFormula is ";

ClosedCurveInterpolation::usage = "ClosedCurveInterpolation[data,t] \
is ";

ShowTimesList::usage = "ShowTimesList is ";

OpenCurveInterpolation::usage = "OpenCurveInterpolation[data,t] is \
";

\(\(PolynomialLevelSetInterpolation::usage = "PolynomialLevelSetInterpolation[data,t,x,y] is ";\)\)

\(\(ClosedCurveInterpolationPlot::usage = "ClosedCurveInterpolationPlot[data] is ";\)\)

\(\(OrientedClosedCurveInterpolationPlot::usage = "OrientedClosedCurveInterpolationPlot[data] is ";\)\)

\(\(FilledClosedCurveInterpolationPlot::usage = "FilledClosedCurveInterpolationPlot[data] is ";\)\)

\(\(OpenCurveInterpolationPlot::usage = "OpenCurveInterpolationPlot[data] is ";\)\)

\(\(OrientedOpenCurveInterpolationPlot::usage = "OrientedOpenCurveInterpolationPlot[data] is ";\)\)

\(\(HamiltonianPlot::usage = "HamiltonianPlot[hamiltonian,{q,qMin,qMax},{p,pMin,pMax},options] draws trajectories of the hamiltonian system with the given hamiltonian. It is based on ContourPlot. Specific options are SymplecticMatrix and ShowField. It accepts also ContourPlot, Graphics`Arrow`, Graphics options.";\)\)

SymplecticMatrix::usage = "SymplecticMatrix is an option for \
HamiltonianPlot specifying what symplectic matrix to use. The default is \
SymplecticMatrix->{{0,1},{-1,0}}.";

NDSolvePlot::usage = "NDSolvePlot[differential equations,options] \
solves the equations with NDSolve and then draws a plot of the solution. For \
a planar system the plot has arrows on it. The output is a couple \
{Graphics,NDSolve output}. It takes all options of NDSolve and \
OrientedParametricPlot.";

\(\(PhasePlot::usage = "PhasePlot[{f[x,y],g[x,y]},{x,xMin,xMax},{y,yMin,yMax},{tMin,tMax},options] makes a phase portrait of the system {x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]]}. Specific options are PlotPoints (that controls how many points there are in the initial values grid), InitialPoints, CurvesWithArrows. It also accepts NDSolve, Graphics`Arrow` and Graphics options.";\)\)

\(\(PoincareMap::usage = 
"PoincareMap[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},options] calculates the Poincare map at time t1 starting at time t0 of the differential system x'=fx, y'=fy on on the region uMin<u<uMax, vMin<v<vMax, where u,v are parameters, and the default parameterization is cartesian (i.e., u=x, v=y). If you want (u,v) to be the polar parameters (r,theta) use the option GridParameterization->Polar. You can specify as GridParameterization any function that takes two variables and gives back a 2-element list, like GridParameterization->Function[{u,v},{u^2,v+u}]. PoincareMap will numerically solve the differential system in a rectangular grid of (u,v) points. The number of points in the two directions is governed by the option GridPoints (that behaves like PlotPoints in ParametricPlot3D), that defaults to 5 in each direction. The Poincaré map in the rest of the region is obtained by ListInterpolation.
Example: PoincareMap[{{x,y},{y,-2x}},{t,0,Pi},{{0,1},{-Pi,Pi}},GridParameterization->Polar].";\)\)

\(\(PoincareMapPlot::usage = 
"PoincareMapPlot[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},options] makes a plot of PoincareMap. The syntax for the non-optional arguments is the same as for PoincareMap. The style of the grid lines is governed by the option GridStyle. Normally you get two plots side-by-side: on the left there is the map at time t0 and on the right the map at time t1. With OnlyEndGrid->True you only get the map at time t1. With StackVertical->True you get the two plots stacked vertically instead of horizontally. With Superimposed->True you get the two plots superimposed.
Example: PoincareMapPlot[{{x,y},{y,-2x}},{t,0,Pi},{{0,1},{-Pi,Pi}},GridParameterization->Polar].";\)\)

\(\(PoincareMapFixedPointDensityPlot::usage = 
"PoincareMapFixedPointDensityPlot[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},epsilon,options] helps in visually searching for fixed points in the PoincareMap. It calculates first the PoincareMap phi, then the function d(x,y)=|phi(x,y)-(x,y)|^2, and makes a DensityPlot of Min[d(x,y),epsilon]. This way the regions where d(x,y) is less than epsilon will be highlighted in grey-black, while the rest of the plot will be white. The plot will be parameterized with (u,v), not with (x,y).
Example: PoincareMapFixedPointDensityPlot[{{x,y},{y,-(x/(1+x^2)+3x(x^2/(1+x^2)))}},{t,0,Pi},{{0,1},{-Pi,Pi}},1/100,GridParameterization->Polar].";\)\)

\(\(PoincareMapFixedPointPlot3D::usage = 
"PoincareMapFixedPointPlot3D[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},{u,u0},{v,v0},options] is similar to PoincareMapFixedPointDensityPlotm, but it makes a ParametricPlot3D instead of a DensityPlot. You will have to adjust PlotRange and BoxRatios to see the minima.
Example: PoincareMapFixedPointPlot3D[{{x,y},{y,-(x/(1+x^2)+3x(x^2/(1+x^2)))}},{t,0,Pi},{{0,1},{-Pi,Pi}},1/100,GridParameterization->Polar].";\)\)

\(\(FindPoincareMapFixedPoint::usage = 
"FindPoincareMapFixedPoint[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},{u,u0},{v,v0},options] applies a FindMinimum to PoincareMap with starting point (u0,v0).
Example: FindPoincareMapFixedPoint[{{x,y},{y,-(x/(1+x^2)+3x(x^2/(1+x^2)))}},{t,0,5},{{.65,.85},{-.1,.1}},{r,.76},{θ,0},GridPoints->{40,11},GridParameterization->Polar].";\)\)

GridPoints::usage = 
"GridPoints is an option for \
PoincareMap, PoincareMapPlot, PoincareMapFixedPointDensityPlot, \
PoincareMapFixedPointPlot3D and FindPoincareMapFixedPoint. It specifies how \
many divisions there will be in the u,v directions, just like PlotPoints does \
for ParametricPlot3D. The default is 5 in each direction.";

\(\(GridParameterization::usage = 
"GridParameterization is an option for PoincareMap, PoincareMapPlot, PoincareMapFixedPointDensityPlot, PoincareMapFixedPointPlot3D and FindPoincareMapFixedPoint. It specifies how the parameters (u,v) are mapped to (x,y). The default is the identity, that is, ({##}&). A predefined value is GridParameterization->Polar, which is shorthand for GridParameterization->Function[{r,theta},r{Cos[theta],Sin[theta]}]. In general as a value for GridParameterizationyou can give any object c that when in the pattern c[x,y], with numeric x,y, gives a list of two numbers {x,y}.";\)\)

Polar::usage = 
"Polar is a possible value for \
the option GridParameterization to be used for PoincareMap, PoincareMapPlot, \
PoincareMapFixedPointDensityPlot, PoincareMapFixedPointPlot3D and \
FindPoincareMapFixedPoint. It is shorthand for \
Function[{r,theta},r{Cos[theta],Sin[theta]}].";

GridStyle::usage = 
"GridStyle is an option for \
PoincareMapPlot and PoincareMapFixedPointDensityPlot. It specifies the style \
of the grid lines. The default is a hue from 0 to 1/2 in one direction and \
from 1/2 to 1 in the other direction. With GridStyle->None all lines will be \
black. To make your own style, just make a list {styleU,styleV} where styleU \
and styleV are each in the format required by ParametricPlot.";

OnlyEndGrid::usage = 
"OnlyEndGrid is an option \
for PoincareMapPlot. With the default OnlyEndGrid->False you will get two \
grids, one at the initial and one at the final time. With OnlyEndGrid->True \
you only get the grid at the final time.";

Superimposed::usage = 
"Superimposed is an option \
for PoincareMapPlot. With the default Superimposed->False you will get two \
grids, one at the initial and one at the final time in separate plots. With \
Superimposed->True you only get the two grid superimposed in one \
plot.";

StackVertical::usage = 
"StackVertical is an \
option for PoincareMapPlot. With the default StackVertical->False you will \
get two grids, one at the initial and one at the final time, side by side. \
With StackVertical->True you get the two grid one below the other.";

SamePlotRange::usage = 
"SamePlotRange is an \
option for PoincareMapPlot. With the default SamePlotRange->True you will get \
two grids, one at the initial and one at the final time, with the same \
PlotRange. With SamePlotRange->False you get Automatic settings for the \
PlotRange in each plot separately.";

Begin["`Private`"];

Arrow stuff

Finding the length of a path:

lineLength[pts_List] := 
        Apply[Plus, 
          Sqrt[Apply[Plus, Rest[pts - RotateRight[pts]]^2, 1]]];

Finding the lengths of the subpaths obtained by taking the first 1,2,3,... points of the given path:

successiveLengths[pts_List] := 
        Map[lineLength, Rest[FoldList[Append, {}, pts]]];

Parameterizing a path by arc length:

byArcLength[pts_List] := 
        With[{witoutDuplicates = 
              pts //. {a___, b_List, b_, c___} :> {a, b, c}}, 
          Function[{t}, 
            Evaluate[
              Through[Map[
                    Interpolation[
                        Transpose[{successiveLengths[witoutDuplicates], #}], 
                        InterpolationOrder -> 1] &, 
                    Transpose[witoutDuplicates]][t]]]]];

Inserting arrow heads along a path in given positions (represented by fractions of the total arc length of the path):

insertArrows[pts_List /; lineLength[pts] > $MachineEpsilon, 
          positions_List, opts___?OptionQ] := 
        Module[{el = lineLength[pts], par = byArcLength[pts]}, {Line[pts], 
            Map[Arrow[par[Max[0, el (# - 1/1000)]], 
                  par[el*Min[# + 1/1000, 1]], FilterOptions[Arrow, opts], 
                  HeadScaling -> Absolute] &, positions]}];

Ignoring arrows when the path is too short:

insertArrows[
          pts_List /; lineLength[pts] ≤  $MachineEpsilon, 
          positions_List, opts___?OptionQ] := {Line[pts]};

Clipping lines and adapting arrows to the display area

Determining whether two points are not in the same (weak) side of the plane with respect to at least one of the straight lines x=xMin, x=xMax, y=yMin, y=yMax:

crossingQ[{xMin_, xMax_}, {yMin_, yMax_}][{{x1_?NumericQ, 
              y1_}, {x2_?NumericQ, y2_}}] := 
        Or[x1 < xMin < x2, x2 < xMin < x1, x1 < xMax < x2, x2 < xMax < x1, 
          y1 < yMin < y2, y2 < yMin < y1, y1 < yMax < y2, 
          y2 < yMax < y1];

When two points make crossingQ true, find the intersection between the earliest found crossed straight line and the segment joining the two points:

addCrossingPoint[{{x1_, y1_}, {x2_, y2_}}, {_, xMax_}, _] /; 
          Or[x1 < xMax < x2, x2 < xMax < x1] := 
        With[{t = -((x2 - xMax)/(x1 - 
                        x2))}, 
{{xMax, 
              t*y1 + (1 - t) y2}, {x2, y2}}];

addCrossingPoint[{{x1_, y1_}, {x2_, y2_}}, {xMin_, _}, _] /; 
          Or[x1 < xMin <  x2, x2 < xMin < x1] := 
        With[{t = -((x2 - xMin)/(x1 - 
                        x2))}, 
{{xMin, 
              t*y1 + (1 - t) y2}, {x2, y2}}];

addCrossingPoint[{{x1_, y1_}, {x2_, y2_}}, _, {_, yMax_}] /; 
          Or[y1 < yMax < y2, y2 < yMax < y1] := 
        With[{t = -((y2 - yMax)/(y1 - 
                        y2))}, 
{{t*x1 + (1 - t) 
                  x2, yMax}, {x2, y2}}];

addCrossingPoint[{{x1_, y1_}, {x2_, y2_}}, _, {yMin_, _}] /; 
          Or[y1 < yMin <  y2, y2 < yMin < y1] := 
        With[{t = -((y2 - yMin)/(y1 - 
                        y2))}, 
{{t*x1 + (1 - t) 
                  x2, yMin}, {x2, y2}}];

addCrossingPoint[cplpts_, lts__] /; 
          Not[crossingQ[lts][cplpts]] := {Last[cplpts]};

Determining where a point is with respect to a rectangle:

strictlyInsideQ[{xMin_, xMax_}, {yMin_, yMax_}][{x_?NumericQ, 
            y_}] := And[xMin < x < xMax, yMin < y < yMax];

weaklyInsideQ[{xMin_, xMax_}, {yMin_, yMax_}][{x_?NumericQ, y_}] := 
        And[xMin ≤  x ≤ xMax, 
          yMin ≤ y ≤ yMax];

weaklyOutsideQ[lts__][pt_List] := 
        Not[strictlyInsideQ[lts][pt]];

strictlyOutsideQ[lts__][pt_List] := 
        Not[weaklyInsideQ[lts][pt]];

Clipping the part of a curve lying outside a rectangle, and then placing arrow heads along it:

clipAndPlaceArrows[gr_, lts__List, inversionCriterion_, 
          arrowDensityParameter_?Positive, opts___?OptionQ] := 
        Module[{strOuts, wklIn, adjClosedPaths, subdivided, clipped, 
            wrongDirectionQ, directedTraj, 
            r}, 
strOuts[a_] := 
            PatternTest[a, 
              Evaluate[strictlyOutsideQ[lts]]]; 
wklIn[a_] := 
            PatternTest[a, 
              Evaluate[
                weaklyInsideQ[
                  lts]]]; 
adjClosedPaths = gr /. 
                  Line[{a_, b___, a_}] :> 
                    provisional[{a, b, a}] /. 
                provisional[{a___, strOuts[b_], c___}] :> 
                  provisional[{b, c, a}] /. 
              provisional -> Line; 
subdivided = 
            adjClosedPaths //. 
              Line[pts_] :> 
                Line[Join[{First[pts]}, 
                    Flatten[
                      Map[addCrossingPoint[#, lts] &, Partition[pts, 2, 1]], 
                      1]]]; 
tst = 
            Or[And[strictlyOutsideQ[lts][#1], strictlyOutsideQ[
                      lts][#2]], 
                And[weaklyInsideQ[lts][#1], weaklyInsideQ[
                      lts][#2]]] &; 
clipped = subdivided /. Line[pts_] :> 
                    Map[Line, 
                      Split[pts, 
                        tst]] /. {Line[{strOuts[__]}] :> 
                    Sequence @@ {}, 
                  Line[{a_}] :> 
                    Sequence @@ {}};; 
directedTraj = 
            clipped /. 
              Line[pts_ /; 
                    And[Length[pts] > 1, 
                      inversionCriterion[pts]]] :> 
                Line[Reverse[pts]]; 
r[pts_] /; 
              lineLength[pts] > $MachineEpsilon := r[pts] = 
              Ceiling[arrowDensityParameter*(lineLength[pts]/
                        Sqrt[Apply[Plus, Apply[Subtract, {lts}, 1]^2]] - 
                      1/15)]; 
r[pts_] /; 
              lineLength[pts] ≤ $MachineEpsilon := 
            1; 
directedTraj /. 
            Line[pts_List] :> 
              insertArrows[pts, 
                Range[1, 2  r[pts] - 1, 2]/(2  Max[r[pts], 1]), 
                FilterOptions[Arrow, opts]]];

OrientedParametricPlot

SetAttributes[OrientedParametricPlot, HoldAll]

Options[OrientedParametricPlot] = 
        Join[{HowManyArrows -> 1, ArrowPositions -> Automatic, 
            ReverseDirectionQ -> False, HeadScaling -> Absolute}, 
          Options[ParametricPlot], 
          DeleteCases[Options[Arrow], HeadScaling -> _]];

OrientedParametricPlot[f_, it_List, opts___?OptionQ] := 
        Module[{howMany, pos, revQ, c, crv, el, par, arrowPos, 
            normalizedArrowPos, 
            arwlst}, 
{howMany, pos, 
              revQ} = {HowManyArrows, ArrowPositions, 
                  ReverseDirectionQ} /. {opts} /. 
              Options[OrientedParametricPlot]; 
c = 
            ParametricPlot[f, it, DisplayFunction -> Identity, 
              Evaluate[
                FilterOptions[ParametricPlot, 
                  opts]]]; 
(arrowPos /; 
                pos =!= Automatic) := 
            Flatten[{pos}]; 
(arrowPos /; 
                pos === Automatic) := 
            Range[1, 2  howMany - 1, 2]/
              Max[2  howMany, 
                1]; 
normalizedArrowPos := (arrowPos - 
                  Min[arrowPos, 0])/(Max[arrowPos, 1] - 
                  Min[arrowPos, 0]); 
Show[
            Graphics[
              c[[1]] /. 
                Line[pts_List] :> 
                  insertArrows[If[revQ, Reverse[pts], pts], 
                    normalizedArrowPos, opts], List @@ Rest[c]], 
            FilterOptions[Graphics, opts], 
            DisplayFunction -> $DisplayFunction]];

Trig Interpolation curves

trigFunctionsAndValueMatrix[conditionsStructure_List, timesList_List, 
          t_Symbol, periodLength_: (2  Pi)] := 
        Module[{numberOfConditions, baseFunctions, f, 
            g}, 
numberOfConditions = conditionsStructure // Flatten // Length; 
baseFunctions[r_] = 
            Take[Flatten[
                Join[{1}, 
                  Table[{Cos[2  Pi*k*r/periodLength], 
                      Sin[2  Pi*k*r/periodLength]}, {k, 
                      Ceiling[numberOfConditions/2] + 1}]]], 
              numberOfConditions + 1]; 
f[
              a_, {k_, j_:  1}] := Derivative[j - 1][g][
              timesList[[k]]]; 
{g[t], 
              Flatten[MapIndexed[f, conditionsStructure, {-1}]]} /. 
            g -> baseFunctions];

Options[
          TrigInterpolation] = {EvaluateOnCoefficients -> Identity, 
          Period -> Automatic};
TrigInterpolation[pointList : {{_?NumericQ, __?NumericQ} .. }, 
          t_Symbol, opts___Rule] := 
        Module[{funct, xlist, derivativesList, xmin, xmax, pl, periodLength, 
            derivativesListStructure, trigValueMatrix, 
            coeff}, 
funct = EvaluateOnCoefficients /. {opts} /. Options[TrigInterpolation]; 
xlist = 
            Map[First, pointList]; 
derivativesList = 
            Map[Rest, pointList]; 
{xmin, 
              xmax} = {Min[xlist], 
              Max[xlist]}; 
pl = Period /. {opts} /. 
              Options[TrigInterpolation]; 
periodLength = 
            Which[pl === Automatic, (xmax - xmin)*3, 
              Element[N[pl], Reals] === True, pl, 
              True, (xmax - 
                    xmin)*3]; 
derivativesListStructure = 
            derivativesList /. _?NumericQ -> 
                1; 
{baseFunctions, trigValueMatrix} = 
            trigFunctionsAndValueMatrix[derivativesListStructure, xlist, t, 
              periodLength]; 
coeff = 
            Map[funct, 
              LinearSolve[trigValueMatrix, 
                Flatten[derivativesList]]]; 
If[
            Head[coeff] === LinearSolve, Null, 
            coeff . baseFunctions]];

Options[TrigInterpolationPlot] = {ShowListPlot -> False, 
          ListPlotStyle -> PointSize[ .02], 
          ShowFormula -> True};
TrigInterpolationPlot[pointList : {{_?NumericQ, __?NumericQ} .. }, 
          xmin_?NumericQ, xmax_?NumericQ, opts___Rule] := 
        Module[{slpQ, sformQ, interpFunct, t, curve, 
            listPl}, 
slpQ = ((ShowListPlot /. {opts} /. Options[TrigInterpolationPlot]) === 
                True); 
sformQ = ((ShowFormula /. {opts} /. Options[TrigInterpolationPlot]) === 
                True); 
interpFunct = 
            TrigInterpolation[pointList, t, 
              FilterOptions[TrigInterpolation, opts]]; 
If[
            interpFunct === Null, Return[Null]]; 
curve = 
            Plot[interpFunct, {t, xmin, xmax}, 
              DisplayFunction -> Identity, 
              Evaluate[
                FilterOptions[Plot, opts]]]; 
listPl := 
            ListPlot[Map[Take[#, 2] &, pointList], 
              PlotStyle -> (ListPlotStyle /. {opts} /. 
                    Options[TrigInterpolationPlot]), 
              DisplayFunction -> Identity]; 
{If[
              slpQ, 
Show[{curve, listPl}, 
                Evaluate[FilterOptions[Graphics, opts]], 
                DisplayFunction -> $DisplayFunction], 
Show[curve, Evaluate[FilterOptions[Graphics, opts]], 
                DisplayFunction -> $DisplayFunction]], 
If[sformQ, 
              Evaluate[interpFunct /. t -> Slot[1]] &]}];

Options[
          ClosedCurveInterpolation] = {EvaluateOnCoefficients -> 
            Identity, 
          ShowTimesList -> False};
ClosedCurveInterpolation[
          pointList : {({_?NumericQ, _?
                        NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. }, 
          t_Symbol, opts___Rule] := 
        Module[{stlQ, funct, x, y, conditionsStructure, conditionsPerPoint, 
            timesList, r, f, g, baseFunctions, trigValueMatrix, coeffx, 
            coeffy, result}, 
stlQ = ((ShowTimesList /. {opts} /. Options[ClosedCurveInterpolation]) === 
                True); 
funct = EvaluateOnCoefficients /. {opts} /. Options[ClosedCurveInterpolation]; 
x = 
            pointList /. {a_?NumericQ, _} :> 
                a; 
y = 
            pointList /. {_, b_?NumericQ} :> 
                b; 
conditionsStructure = 
            x /. _?NumericQ -> 
                1; 
conditionsPerPoint = 
            Map[Length[Flatten[#]] &, pointList]/
              2; 
timesList = 
            2  Pi (Drop[#, -1]/Last[#] &)[
                FoldList[Plus, 0, 
                  ListConvolve[{1, 1}, 
                    conditionsPerPoint, {1, 
                      1}]]]; 
{baseFunctions, 
              trigValueMatrix} = 
            trigFunctionsAndValueMatrix[conditionsStructure, timesList, 
              t]; 
coeffx = 
            Map[funct, 
              LinearSolve[trigValueMatrix, 
                Flatten[x]]]; 
coeffy = 
            Map[funct, 
              LinearSolve[trigValueMatrix, 
                Flatten[
                  y]]]; 
result = 
{baseFunctions . coeffx, 
              baseFunctions . coeffy}; 
If[
            Or[Head[coeffx] === LinearSolve, Head[coeffy] === LinearSolve], 
            Null, {result, 
              If[stlQ, timesList, Sequence @@ {}]}]];

Options[
          OpenCurveInterpolation] = {EvaluateOnCoefficients -> Identity, 
          ShowTimesList -> False, 
          Period -> Automatic};
OpenCurveInterpolation[
          pointList : {({_?NumericQ, _?
                        NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. }, 
          t_Symbol, opts___Rule] := 
        Module[{stlQ, funct, x, y, conditionsStructure, conditionsPerPoint, 
            timesList, tmax, pl, periodLength, baseFunctions, 
            trigValueMatrix, coeffx, coeffy, 
            result}, 
stlQ = ((ShowTimesList /. {opts} /. Options[ClosedCurveInterpolation]) === 
                True); 
funct = EvaluateOnCoefficients /. {opts} /. Options[ClosedCurveInterpolation]; 
x = 
            pointList /. {a_?NumericQ, _} :> 
                a; 
y = 
            pointList /. {_, b_?NumericQ} :> 
                b; 
conditionsStructure = 
            x /. _?NumericQ -> 
                1; 
conditionsPerPoint = 
            Map[Length[Flatten[#]] &, pointList]/
              2; 
timesList = 
            FoldList[Plus, 0, ListConvolve[{1, 1}, conditionsPerPoint]]; 
          tmax = Last[
              timesList]; 
pl = Period /. {opts} /. 
              Options[OpenCurveInterpolation]; 
periodLength = 
            Which[pl === Automatic, tmax*3, Element[N[pl], Reals] === True, 
              pl, True, tmax*3]; 
{baseFunctions, 
              trigValueMatrix} = 
            trigFunctionsAndValueMatrix[conditionsStructure, timesList, t, 
              periodLength]; 
coeffx = 
            Map[funct, 
              LinearSolve[trigValueMatrix, 
                Flatten[x]]]; 
coeffy = 
            Map[funct, 
              LinearSolve[trigValueMatrix, 
                Flatten[
                  y]]]; 
result = 
{baseFunctions . coeffx, 
              baseFunctions . coeffy}; 
If[
            Or[Head[coeffx] === LinearSolve, Head[coeffy] === LinearSolve], 
            Null, {result, 
              If[stlQ, timesList, Sequence @@ {}]}]];

Options[PolynomialLevelSetInterpolation] = 
        Options[ClosedCurveInterpolation];

\(\(PolynomialLevelSetInterpolation::"NonExactInput" = "PolynomialLevelSetInterpolation is designed only for exact (infinite precision) numeric input.";\)\)

PolynomialLevelSetInterpolation[
          pointList : {({_?NumericQ, _?
                          NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. } /; 
            Precision[pointList] < Infinity, t_Symbol, x_Symbol, y_Symbol, 
          opts___Rule] := (Message[
            PolynomialLevelSetInterpolation::"NonExactInput"]; 
          Null);

PolynomialLevelSetInterpolation[
          pointList : {({_?NumericQ, _?
                          NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. } /; 
            Precision[pointList] == Infinity, t_Symbol, x_Symbol, 
          y_Symbol, opts___Rule] := 
        Module[{funct, parametricContour, poly, a, b, 
            simplifiedPoly}, 
funct = EvaluateOnCoefficients /. {opts} /. 
              Options[ClosedCurveInterpolation]; 
parametricContour = 
            ClosedCurveInterpolation[pointList, t, 
              FilterOptions[ClosedCurveInterpolation, 
                opts]]; 
If[parametricContour === Null, 
            Return[Null]]; 
          poly = Eliminate[
                Join[Thread[{x, y} == 
                        TrigExpand[
                          Take[Flatten[parametricContour], 2]] /. {Cos[
                            t] -> a, Sin[t] -> b}], {a^2 + 
                        b^2 == 1}], {a, b}] /. 
              Equal -> Subtract; 
simplifiedPoly = 
            Collect[poly, {x^(_ : 1) y^(_ : 1), x^(_ : 1), 
                y^(_ : 1)}, 
              funct]; 
{parametricContour, 
            simplifiedPoly}];

Options[ClosedCurveInterpolationPlot] = 
      Join[{ShowListPlot -> False, 
          ListPlotStyle -> PointSize[ .02], ShowFormula -> True}, 
        Options[ParametricPlot]]; 
    ClosedCurveInterpolationPlot[
        pointList : {({_?NumericQ, _?
                      NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. }, 
        opts___] := 
      Module[{slpQ, sformQ, t, paramCurve, basicPlot, 
          pointsPlot}, 
slpQ = ((ShowListPlot /. {opts} /. Options[ClosedCurveInterpolationPlot]) === 
              True); 
sformQ = ((ShowFormula /. {opts} /. Options[ClosedCurveInterpolationPlot]) === 
              True); 
paramCurve = 
          ClosedCurveInterpolation[pointList, t, 
            Evaluate[
              FilterOptions[ClosedCurveInterpolation, 
                opts]]]; 
If[paramCurve === Null, 
          Return[Null]]; 
basicPlot = 
          ParametricPlot[Evaluate[paramCurve[[1]]], {t, 0, 2  Pi}, 
            DisplayFunction -> Identity, 
            Evaluate[
              FilterOptions[ParametricPlot, 
                opts]]]; 
pointsPlot := 
          ListPlot[Map[Take[Flatten[#], 2] &, pointList], 
            DisplayFunction -> Identity, 
            PlotStyle -> (ListPlotStyle /. {opts} /. 
                  Options[
                    ClosedCurveInterpolationPlot])]; 
{Show[
            If[slpQ, {basicPlot, 
                pointsPlot}, 
basicPlot], 
Evaluate[FilterOptions[Graphics, opts]], 
            DisplayFunction -> $DisplayFunction], 
          If[sformQ, {Evaluate[paramCurve[[1]] /. t -> Slot[1]] &, 
              Sequence @@ Rest[paramCurve]}], Sequence @@ {}}];

Options[OrientedClosedCurveInterpolationPlot] = 
      Join[{ShowListPlot -> False, 
          ListPlotStyle -> PointSize[ .02], ShowFormula -> True, 
          HowManyArrows -> 3}, Options[OrientedParametricPlot]]; 
    OrientedClosedCurveInterpolationPlot[
        pointList : {({_?NumericQ, _?
                      NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. }, 
        opts___] := 
      Module[{slpQ, sformQ, t, paramCurve, basicPlot, 
          pointsPlot}, 
slpQ = ((ShowListPlot /. {opts} /. Options[ClosedCurveInterpolationPlot]) === 
              True); 
sformQ = ((ShowFormula /. {opts} /. Options[ClosedCurveInterpolationPlot]) === 
              True); 
paramCurve = 
          ClosedCurveInterpolation[pointList, t, 
            Evaluate[
              FilterOptions[ClosedCurveInterpolation, 
                opts]]]; 
If[paramCurve === Null, 
          Return[Null]]; 
basicPlot = 
          OrientedParametricPlot[
            Evaluate[paramCurve[[1]]], {t, 0, 2  Pi}, 
            DisplayFunction -> Identity, 
            Evaluate[
              FilterOptions[OrientedParametricPlot, opts, 
                Sequence @@ Options[OrientedClosedCurveInterpolationPlot]]], 
            Evaluate[
              FilterOptions[ParametricPlot, 
                opts]]]; 
pointsPlot := 
          ListPlot[Map[Take[Flatten[#], 2] &, pointList], 
            DisplayFunction -> Identity, 
            PlotStyle -> (ListPlotStyle /. {opts} /. 
                  Options[
                    ClosedCurveInterpolationPlot])]; 
{Show[
            If[slpQ, {basicPlot, 
                pointsPlot}, 
basicPlot], 
Evaluate[FilterOptions[Graphics, opts]], 
            DisplayFunction -> $DisplayFunction], 
          If[sformQ, {Evaluate[paramCurve[[1]] /. t -> Slot[1]] &, 
              Sequence @@ Rest[paramCurve]}], Sequence @@ {}}];

Options[
        FilledClosedCurveInterpolationPlot] = {ShowListPlot -> False, 
        ListPlotStyle -> PointSize[ .02], ShowFormula -> True, 
        Fills -> GrayLevel[ .5], Curves -> Front}; 
    FilledClosedCurveInterpolationPlot[
        pointList : {({_?NumericQ, _?
                      NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. }, 
        opts___] := 
      Module[{slpQ, sformQ, crv, crvsQ, fls, t, basicPlot, paramCurve, 
          filledPolygon, 
          pointsPlot}, 
slpQ = ((ShowListPlot /. {opts} /. Options[ClosedCurveInterpolationPlot]) === 
              True); 
sformQ = ((ShowFormula /. {opts} /. Options[ClosedCurveInterpolationPlot]) === 
              True); 
crv = (Curves /. {opts} /. 
              Options[ClosedCurveInterpolationPlot]); 
crvsQ = 
          Which[MemberQ[{True, Front}, crv], True, 
            MemberQ[{False, None}, crv], False, True, 
            True]; 
fls = Fills /. {opts} /. 
            Options[FilledClosedCurveInterpolationPlot]; 
crv = 
          ClosedCurveInterpolationPlot[pointList, 
            DisplayFunction -> Identity, 
            Evaluate[FilterOptions[ClosedCurveInterpolation, opts]], 
            Evaluate[FilterOptions[ClosedCurveInterpolationPlot, opts]], 
            Evaluate[
              FilterOptions[ParametricPlot, opts]]]; 
If[
          crv === Null, 
          Return[Null]]; 
{basicPlot, paramCurve} = 
          crv; 
filledPolygon = 
          Graphics[
            Cases[basicPlot, 
              Line[pts_List] :> {fls, Polygon[pts]}, {1, 
                Infinity}]]; 
pointsPlot := 
          ListPlot[Map[Take[Flatten[#], 2] &, pointList], 
            DisplayFunction -> Identity, 
            PlotStyle -> (ListPlotStyle /. {opts} /. 
                  Options[
                    ClosedCurveInterpolationPlot])]; 
{Show[{filledPolygon, 
              If[crvsQ, basicPlot, Sequence @@ {}], 
              If[slpQ, pointsPlot, 
                Sequence @@ {}]}, 
Evaluate[
              FilterOptions[Graphics, opts]], 
            DisplayFunction -> $DisplayFunction], 
          If[sformQ, {Evaluate[paramCurve[[1]]], 
              Sequence @@ Rest[paramCurve]}], Sequence @@ {}}];

Options[OpenCurveInterpolationPlot] = 
      Join[{ShowListPlot -> False, 
          ListPlotStyle -> PointSize[ .02], ShowFormula -> True}, 
        Options[ParametricPlot]]; 
    OpenCurveInterpolationPlot[
        pointList : {({_?NumericQ, _?
                      NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. }, 
        opts___] := 
      Module[{slpQ, sformQ, stlstQ, t, paramCurve, tmax, basicPlot, 
          pointsPlot}, 
slpQ = ((ShowListPlot /. {opts} /. Options[OpenCurveInterpolationPlot]) === 
              True); 
sformQ = ((ShowFormula /. {opts} /. Options[OpenCurveInterpolationPlot]) === 
              True); 
stlstQ = ((ShowTimesList /. {opts} /. Options[OpenCurveInterpolationPlot]) === 
              True); 
paramCurve = 
          OpenCurveInterpolation[pointList, t, ShowTimesList -> True, 
            Evaluate[
              FilterOptions[OpenCurveInterpolation, 
                opts]]]; 
If[paramCurve === Null, 
          Return[Null]]; 
tmax = 
          paramCurve[[-1, -1]]; 
paramCurve = 
          If[stlstQ === True, 
            paramCurve, {First[paramCurve]}]; 
If[
          paramCurve === Null, Return[Null]]; 
basicPlot = 
          ParametricPlot[Evaluate[paramCurve[[1]]], {t, 0, tmax}, 
            DisplayFunction -> Identity, 
            Evaluate[
              FilterOptions[ParametricPlot, 
                opts]]]; 
pointsPlot := 
          ListPlot[Map[Take[Flatten[#], 2] &, pointList], 
            DisplayFunction -> Identity, 
            PlotStyle -> (ListPlotStyle /. {opts} /. 
                  Options[
                    OpenCurveInterpolationPlot])]; 
{Show[
            If[slpQ, {basicPlot, 
                pointsPlot}, 
basicPlot], 
Evaluate[FilterOptions[Graphics, opts]], 
            DisplayFunction -> $DisplayFunction], 
          If[sformQ, {Evaluate[paramCurve[[1]] /. t -> Slot[1]] &, 
              Sequence @@ Rest[paramCurve]}], Sequence @@ {}}];

Options[OrientedOpenCurveInterpolationPlot] = 
      Join[{ShowListPlot -> False, 
          ListPlotStyle -> PointSize[ .02], ShowFormula -> True}, 
        Options[OrientedParametricPlot]]; 
    OrientedOpenCurveInterpolationPlot[
        pointList : {({_?NumericQ, _?
                      NumericQ} | {{_?NumericQ, _?NumericQ} .. }) .. }, 
        opts___] := 
      Module[{slpQ, sformQ, stlstQ, t, paramCurve, tmax, basicPlot, 
          pointsPlot}, 
slpQ = ((ShowListPlot /. {opts} /. Options[OpenCurveInterpolationPlot]) === 
              True); 
sformQ = ((ShowFormula /. {opts} /. Options[OpenCurveInterpolationPlot]) === 
              True); 
stlstQ = ((ShowTimesList /. {opts} /. Options[OpenCurveInterpolationPlot]) === 
              True); 
paramCurve = 
          OpenCurveInterpolation[pointList, t, ShowTimesList -> True, 
            Evaluate[
              FilterOptions[OpenCurveInterpolation, 
                opts]]]; 
If[paramCurve === Null, 
          Return[Null]]; 
tmax = 
          paramCurve[[-1, -1]]; 
paramCurve = 
          If[stlstQ === True, 
            paramCurve, {First[paramCurve]}]; 
basicPlot = 
          OrientedParametricPlot[Evaluate[paramCurve[[1]]], {t, 0, tmax}, 
            DisplayFunction -> Identity, 
            Evaluate[
              FilterOptions[OrientedParametricPlot, 
                opts]]]; 
pointsPlot := 
          ListPlot[Map[Take[Flatten[#], 2] &, pointList], 
            DisplayFunction -> Identity, 
            PlotStyle -> (ListPlotStyle /. {opts} /. 
                  Options[
                    OpenCurveInterpolationPlot])]; 
{Show[
            If[slpQ, {basicPlot, 
                pointsPlot}, 
basicPlot], 
Evaluate[FilterOptions[Graphics, opts]], 
            DisplayFunction -> $DisplayFunction], 
          If[sformQ, {Evaluate[paramCurve[[1]] /. t -> Slot[1]] &, 
              Sequence @@ Rest[paramCurve]}], Sequence @@ {}}];

HamiltonianPlot

Options[
          HamiltonianPlot] = {SymplecticMatrix -> {{0, 1}, {-1, 0}}, 
          ShowField -> False, ArrowDensity -> 1};

HamiltonianPlot[
          h_, {x_Symbol, xMin_, xMax_}, {y_Symbol, yMin_, yMax_}, 
          opts___?OptionQ] := 
        Module[{levelData, sympMat, showFieldQ, arwDnsty, approx, grad, a, b, 
            wrongDirectionQ, allTrajs, lateralMargin, topBottomMargin, 
            fld}, 
levelData = 
            ContourPlot[h, {x, xMin, xMax}, {y, yMin, yMax}, 
              DisplayFunction -> Identity, 
              Evaluate[FilterOptions[ContourPlot, opts]], 
              ContourShading -> False]; 
{sympMat, 
              showFieldQ, 
              arwDnsty} = {SymplecticMatrix, ShowField, 
                  ArrowDensity} /. {opts} /. 
              Options[HamiltonianPlot]; 
approx = 
            ListInterpolation[
              Transpose[
                levelData[[1]]], {{xMin, xMax}, {yMin, 
                  yMax}}]; 
grad[{a_, b_}] := 
            Through[{Derivative[1, 0][approx], Derivative[0, 1][
                    approx]}[a, b]]; 
wrongDirectionQ[
              pts_List] := (byArcLength[pts]'[
                    0] . (sympMat . grad[First[pts]]) < 
                0); 
allTrajs = 
            clipAndPlaceArrows[
              Graphics[levelData], {xMin, xMax}, {yMin, yMax}, 
              wrongDirectionQ, arwDnsty, 
              opts]; 
{lateralMargin, 
              topBottomMargin} = {xMax - xMin, yMax - yMin}/
              40; 
fld := 
            PlotVectorField[
              sympMat . grad[{a, b}], {a, xMin + lateralMargin, 
                xMax - lateralMargin}, {b, yMin + topBottomMargin, 
                yMax - topBottomMargin}, 
              DisplayFunction -> 
                Identity]; 
Show[{allTrajs, 
              If[showFieldQ, fld, Sequence @@ {}]}, 
            FilterOptions[Graphics, opts], 
            DisplayFunction -> $DisplayFunction]];

NDSolvePlot

The case of a system with two unknown functions:

NDSolvePlot[eqs_, {x1_, x2_}, {t_, tMin_, tMax_}, opts___?OptionQ] := 
        Module[{solRule, sol, f}, 
solRule = 
            NDSolve[eqs, {x1, x2}, {t, tMin, tMax}, 
              FilterOptions[NDSolve, opts]]; 
sol = {x1, 
                x2} /. solRule[[1]]; 
f = 
            If[Head[sol[[1]]] === InterpolatingFunction, Through[sol[t]], 
              sol]; 
{solRule, 
            OrientedParametricPlot[Evaluate[f], {t, tMin, tMax}, 
              Evaluate[
                FilterOptions[OrientedParametricPlot, opts]]]}];

The case of a single unknown function:

NDSolvePlot[eqs_, x_ /; Head[x] =!= List, {t_, tMin_, tMax_}, 
          opts___?OptionQ] := 
        Module[{solRule, sol, f}, 
          solRule = 
            NDSolve[eqs, x, {t, tMin, tMax}, 
              FilterOptions[NDSolve, opts]]; 
sol = 
            x /. solRule[[1]]; 
f = 
            If[Head[sol] === InterpolatingFunction, sol[t], 
              sol]; 
{solRule, 
            OrientedParametricPlot[Evaluate[{t, f}], {t, tMin, tMax}, 
              Evaluate[FilterOptions[OrientedParametricPlot, opts]], 
              Evaluate[FilterOptions[Plot, opts]], 
              HowManyArrows -> 0]}];

PhasePlot

Options[PhasePlot] = {PlotPoints -> 5, 
          InitialPoints -> Automatic, CurvesWithArrows -> True, 
          Frame -> True,  ArrowDensity -> 1, 
          ShowField -> False};

SetAttributes[PhasePlot, HoldAll];

PhasePlot::"BadPlotPoints" = "Bad PlotPoints";

PhasePlot[
          field_List, {x_Symbol, xMin_, xMax_} /; 
            xMin < xMax, {y_Symbol, yMin_, yMax_} /; 
            yMin < yMax, {tMin_?NonPositive, tMax_?NonNegative}, 
          opts___?OptionQ] := 
        Module[{hm, shwFldQ, shwArrsQ, arwDnsty, lateralMargin, 
            topBottomMargin, fld, ptsX, ptsY, intlPts, strtPts, eq, t, u, v, 
            x0, y0, sol, allSols, 
            allTrajs}, 
{hm, shwFldQ, shwArrsQ, 
              arwDnsty} = {PlotPoints, ShowField, CurvesWithArrows, 
                  ArrowDensity} /. {opts} /. 
              Options[PhasePlot]; 
{lateralMargin, 
              topBottomMargin} = {xMax - xMin, yMax - yMin}/
              40; 
fld := 
            PlotVectorField[
              field, {x, xMin + lateralMargin, xMax - lateralMargin}, {y, 
                yMin + topBottomMargin, yMax - topBottomMargin}, 
              DisplayFunction -> Identity]; 
{ptsX, 
              ptsY} = Which[
              MatchQ[hm, {_Integer?Positive, __Integer?Positive}], hm, 
              MatchQ[hm, _Integer?Positive], {hm, hm}, True, 
              Message[PhasePlot::BadPlotPoints]; 
              Abort[]]; 
intlPts = InitialPoints /. {opts} /. Options[PhasePlot]; 
strtPts = 
            If[intlPts === Automatic, 
              Flatten[Outer[List, 
                  xMin + (xMax - xMin) 
                      Range[1, 2  ptsX - 1, 2]/(2  ptsX), 
                  yMin + (yMax - yMin) 
                      Range[1, 2  ptsY - 1, 2]/(2  ptsY)], 1], 
              intlPts]; 
eq[x0_, 
              y0_] := {{u'[t] == field[[1]], v'[t] == 
                    field[[2]], u[0] == x0, 
                  v[0] == y0} /. {x -> u[t], y -> v[t]}, {u, 
                v}, {t, tMin, tMax}, DisplayFunction -> Identity, 
              HowManyArrows -> 0, 
              FilterOptions[NDSolve, opts]}; 
sol[{x0_, 
                y0_}] := Apply[NDSolvePlot, eq[x0, y0]][[2, 1, 1, 1, 
                1]]; 
allSols = 
            Map[sol, strtPts]; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[allSols, {xMin, xMax}, {yMin, yMax}, 
                False &, arwDnsty, opts], 
              allSols]; 
Show[{If[shwFldQ, fld, 
                Sequence @@ {}], Graphics[allTrajs]}, 
            FilterOptions[Graphics, opts], 
            PlotRange -> {{xMin, xMax}, {yMin, yMax}}, 
            FilterOptions[Graphics, 
              Sequence @@ Options[PhasePlot]]]];

LinearPhasePlot basic tests

Error messages when a bad matrix is given:

admissibleMatrixQ::"notMatrix" = "Not a 2 by 2 matrix";

admissibleMatrixQ::"nonNumeric" = "Not a numeric \
matrix";

admissibleMatrixQ::"complex" = "Complex matrix";

admissibleMatrixQ::"inexactMatrix" = "Not an exact numeric \
matrix";

Determining whether the given matrix is good:

admissibleMatrixQ[mat_?MatrixQ] := admissibleMatrixQ[mat] = 
          And[
If[
              MatrixQ[mat] && Length[mat] === 2 && TensorRank[mat] === 2, 
              True, Message[admissibleMatrixQ::"notMatrix"]; 
              False], 
If[
              And @@ Map[NumericQ, Flatten[mat]], True, 
              Message[admissibleMatrixQ::"nonNumeric"]; 
              False], 
If[
              And @@ Map[FreeQ[Chop[N[#]], Complex] &, Flatten[mat]], True, 
              Message[admissibleMatrixQ::"complex"]; False], 
            If[Precision[mat] === Infinity, True, 
              Message[admissibleMatrixQ::"inexactMatrix"]; 
              False]];

The simplest matrix structure tests:

diagonalQ[mat_?MatrixQ] := 
        Simplify[mat[[1, 2]]] === 0 && 
          Simplify[mat[[2, 1]]] === 0;

nonNullQ[mat_] := 
        Apply[Or, Map[# =!= 0 &, Simplify[Flatten[{mat}]]]];

nullMatrixQ[mat_] := 
        Union[Flatten[Simplify[mat]]] === {0};

LinearPhasePlot basic matrix functions

The trace of a matrix requires to test what Mathematica version is running:

matrixTrace = 
        If[$VersionNumber < 4, (#[[1, 1]] + #[[2, 2]]) &, 
          Tr];

The discriminant of a 2 by 2 real matrix:

discriminant[mat_] := 
        Simplify[matrixTrace[mat]^2 - 4  Det[mat]];

The real canonical form of a real 2 by 2 matrix:

realCanonicalForm[mat_ /; discriminant[mat] ≥ 0] := 
        JordanDecomposition[m];

realCanonicalForm[mat_ /; discriminant[mat] < 0] := 
        Module[{complexSim, complexJordan, 
            realSim}, 
{complexSim, complexJordan} = 
            JordanDecomposition[mat]; 
realSim = 
            Transpose[{Re[#], Im[#]} &[complexSim // Transpose // 
                  First]]; 
{realSim, 
            Inverse[realSim] . mat . realSim}];

LinearPhasePlot spectral analysis

The spectral classification of 2 by 2 matrices:

saddleQ[mat_?admissibleMatrixQ] := 
        discriminant[mat] > 0 && Simplify[Det[mat]] < 0;

nodeQ[mat_?admissibleMatrixQ] := 
        discriminant[mat] > 0 && Simplify[Det[mat]] > 0;

singularDistinctEigenvaluesQ[mat_?admissibleMatrixQ] := 
        discriminant[mat] > 0 && Simplify[Det[mat]] === 0;

improperNodeQ[mat_?admissibleMatrixQ] := 
        discriminant[mat] === 0 && 
          Simplify[matrixTrace[mat]] =!= 
            0 && ! diagonalQ[mat];

nonNullMultipleOfIdentity[mat_?admissibleMatrixQ] := 
        discriminant[mat] === 0 && diagonalQ[mat] && 
          nonNullQ[mat];

nonNullNilpotentQ[mat_?admissibleMatrixQ] := 
        discriminant[mat] === 0 && Simplify[matrixTrace[mat]] === 0 && 
          nonNullQ[mat];

centerQ[mat_?admissibleMatrixQ] := 
        discriminant[mat] < 0 && 
          Simplify[matrixTrace[mat]] === 0;

spiralQ[mat_?admissibleMatrixQ] := 
        discriminant[mat] < 0 && 
          Simplify[matrixTrace[mat]] =!= 0;

LinearPhasePlot main

preamble

Options[
          LinearPhasePlot] = {PlotRange -> {{-1, 1}, {-1, 1}}, 
          AspectRatio -> Automatic, Frame -> True, 
          FrameTicks -> {{0}, {0}}, 
          NotableLineStyle -> RGBColor[1, 0, 0], 
          NumberOfTrajectories -> 5, WriteReport -> True, 
          ShowField -> False, CurvesWithArrows -> True,  
          ArrowDensity -> 1};

inverseRange[s_?MatrixQ] := 
        Map[Max, 
          Inverse[s] . {{1, -1, -1, 1}, {1, 
                1, -1, -1}}];

equil = {PointSize[ .03], Point[{0, 0}], GrayLevel[1], 
          PointSize[ .02], Point[{0, 0}]};

wrongDirectionQ[mat_?MatrixQ][{pt1_List, 
            pt2_, ___}] := ((pt2 - pt1) . (mat . (pt1 + pt2)/
                  2) < 0);

Node (two nonvanishing real eigenvalues of the same sign)

LinearPhasePlot[mat_?nodeQ, opts___Rule] := 
        Module[{sim, diag, λ1, λ2, wrtRptQ, shwFldQ, 
            shwArrsQ, arwDnsty, n, style, fld, xmax, ymax, sol, vert, orizz, 
            r, coeffList, basicTrajs, 
            allTrajs}, 
{sim, diag} = 
            JordanDecomposition[
              mat]; 
{λ1, λ2} = 
            Transpose[diag, {1, 1}]; 
{wrtRptQ, shwFldQ, 
              shwArrsQ, arwDnsty, n, 
              style} = {WriteReport, ShowField, CurvesWithArrows,  
                  ArrowDensity, NumberOfTrajectories + 1, 
                  NotableLineStyle} /. {opts} /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print[If[λ1 > 
                    0, "Repulsive", "Attractive"] <> " node"], 
            Null]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, - .9,  .9}, {y, - .9,  .9}, 
              DisplayFunction -> Identity]; 
{xmax, 
              ymax} = N[inverseRange[sim]]; 
sol[c_] := 
            Map[N[sim] . # &, Plot[c*Abs[x]^(λ2/λ1), 
                    Evaluate[{x, -Min[xmax, r], Min[xmax, r]} /. 
                        r :> (ymax/
                                Abs[c])^(λ1/λ2)], 
                    DisplayFunction -> Identity][[1, 1, 1]] /. 
                Line[{a___, {x1_?NonPositive, y1_}, {x2_?NonNegative, y2_}, 
                      b___}] :> {Line[{a, {x1, y1}, {0, 0}}], 
                    Line[{{0, 0}, {x2, y2}, 
                        b}]}, {3}]; 
vert = {style, 
                Map[N[sim] . # &, {Line[{{0, -ymax}, -{0, 0}/10^3}], 
                    Line[{{0, 0}/10^3, {0, ymax}}]}, {3}]} // 
              Flatten; 
orizz = {style, 
                Map[N[sim] . # &, {Line[{{-xmax, 0}, { .01,  .01}}], 
                    Line[{{0, 0}/10^3, {xmax, 0}}]}, {3}]} // 
              Flatten; 
coeffList = Join[#, -#] &[
              Table[(1 - r)*
                  ymax/(r*xmax)^(λ2/λ1), {r, 1/n, 
                  1 - 1/n, 1/n}]]; 
basicTrajs = {Map[sol, 
                coeffList], vert, orizz, 
              equil}; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[basicTrajs, {-1, 1}, {-1, 1}, 
                wrongDirectionQ[mat], arwDnsty, opts], 
              basicTrajs]; 
Show[{If[shwFldQ, fld, 
                Sequence @@ {}], Graphics[allTrajs]}, 
            PlotRange -> {{-1, 1}, {-1, 1}}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];

Saddle (two nonvanishing real eigenvalues with opposite signs)

\(\(LinearPhasePlot[mat_?saddleQ, opts___Rule] := 
        Module[{sim, diag, λ1, λ2, wrtRptQ, shwFldQ, 
            shwArrsQ, arwDnsty, n, style, fld, xmax, ymax, sol0, sol, 
            coeffList, vert, orizz, basicTrajs, 
            allTrajs}, 
{sim, diag} = 
            JordanDecomposition[
              mat]; 
{λ1, λ2} = {diag[\([1, 
                  1]\)], diag[\([2, 2]\)]}; 
{wrtRptQ, 
              shwFldQ, shwArrsQ, arwDnsty, n, 
              style} = \({WriteReport, ShowField, CurvesWithArrows, \ 
                  ArrowDensity, NumberOfTrajectories + 1, 
                  NotableLineStyle} /. {opts}\) /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print["Saddle"], Null]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, \(- .9\),  .9}, {y, \(- .9\),  .9}, 
              DisplayFunction -> Identity]; 
{xmax, 
              ymax} = N[inverseRange[sim]]; 
sol0[
              c_] := \(Plot[c*Abs[x]^\((λ2/λ1)\), 
                Evaluate[{x, Abs[ymax/c]^\((λ1/λ2)\), 
                    Min[xmax, 1]}], DisplayFunction -> Identity]\)[\([1, 
                1, 1]\)]; 
sol[
              c_] := {Map[N[sim] . # &, sol0[c], {2}], 
              Map[N[sim] . {{\(-1\), 0}, {0, 1}} . # &, 
                sol0[c], {2}]}; 
coeffList = \(Join[#, \(-#\)] &\)[Table[r*ymax/\((r*xmax)\)^\((λ2/λ1)\), {r, 1/n, 
                  1 - 1/n, 1/n}]]; 
vert = {style, 
                Map[N[sim] . # &, {Line[{{0, \(-ymax\)}, {0, 0}}], 
                    Line[{{0, 0}, {0, ymax}}]}, {3}]} // 
              Flatten; 
orizz = {style, 
                Map[N[sim] . # &, {Line[{{\(-xmax\), 0}, {0, 0}}], 
                    Line[{{0, 0}, {xmax, 0}}]}, {3}]} // 
              Flatten; 
basicTrajs = {Map[sol, coeffList], 
              vert, orizz, equil}; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[basicTrajs, {\(-1\), 1}, {\(-1\), 1}, 
                wrongDirectionQ[mat], arwDnsty, opts], basicTrajs]; 
          Show[{If[shwFldQ, fld, Sequence @@ {}], Graphics[allTrajs]}, 
            PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];\)\)

Two distinct real eigenvalues, one vanishing

LinearPhasePlot[mat_?singularDistinctEigenvaluesQ, opts___?OptionQ] := 
        Module[{sim, canon, λ, wrtRptQ, shwFldQ, shwArrsQ, arwDnsty, 
            n, style, fld, trasf, xmax, ymax, sol, vert, coeffList, 
            basicTrajs, allTrajs}, 
{sim, canon} = 
            JordanDecomposition[
              mat]; 
λ = Select[
                Flatten[canon], # =!= 0 &, 
                1][[1]]; 
{wrtRptQ, shwFldQ, 
              shwArrsQ, arwDnsty, n, 
              style} = {WriteReport, ShowField, CurvesWithArrows,  
                  ArrowDensity, NumberOfTrajectories, 
                  NotableLineStyle} /. {opts} /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print["One null eigenvalue, the other one " <> 
                If[λ > 0, "positive", "negative"]], 
            Null]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, - .9,  .9}, {y, - .9,  .9}, 
              DisplayFunction -> Identity]; 
trasf = 
            N[sim] . 
                If[canon[[1, 1]] === λ, 
                  IdentityMatrix[
                    2], {{0, 1}, {1, 0}}] . # &; 
{xmax, 
              ymax} = N[inverseRange[sim]]; 
sol[c_] := 
            Map[trasf, {Line[{{-xmax, c}, {0, c}}], 
                Line[{{0, c}, {xmax, 
                      c}}]}, {3}]; 
vert = {{GrayLevel[1], 
                Thickness[ .04], 
                Line[trasf /@ {{0, -ymax}, {0, ymax}}]}, {RGBColor[1, 0, 
                    0], Dashing[{ .02,  .02}], Thickness[ .01], style, 
                  Line[trasf /@ {{0, -ymax}, {0, ymax}}]} // 
                Flatten}; 
coeffList = 
            ymax*Range[-n, n]/n; 
basicTrajs = 
            Map[sol, coeffList]; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[basicTrajs, {-1, 1}, {-1, 1}, 
                wrongDirectionQ[mat], arwDnsty, opts], basicTrajs]; 
          Show[{If[shwFldQ, fld, Sequence @@ {}], 
              Graphics[{allTrajs, vert, equil}]}, 
            PlotRange -> {{-1, 1}, {-1, 1}}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];

Improper node (two real coincident eigenvalues, but the matrix is not diagonal)

LinearPhasePlot[mat_?improperNodeQ, opts___Rule] := 
        Module[{sim, canon, λ, wrtRptQ, shwFldQ, shwArrsQ, arwDnsty, 
            n, style, fld, xmax, ymax, sol, u, coeffList, orizz, basicTrajs, 
            allTrajs}, 
{sim, canon} = 
            N[JordanDecomposition[mat]]; 
λ = 
            canon[[1, 1]]; 
{wrtRptQ, shwFldQ, 
              shwArrsQ, arwDnsty, n, 
              style} = {WriteReport, ShowField, CurvesWithArrows,  
                  ArrowDensity, 2  NumberOfTrajectories, 
                  NotableLineStyle} /. {opts} /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print[If[λ > 
                    0, "Repulsive", "Attractive"] <> " improper \
node"], Null]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, - .9,  .9}, {y, - .9,  .9}, 
              DisplayFunction -> Identity]; 
{xmax, 
              ymax} = N[inverseRange[sim]]; 
sol[c_] := 
            Map[N[sim] . 
                  Reverse[#] &, Plot[
                    c*y + y*
                        Log[Max[
                              Abs[y],  .000001  
                                ymax]]/λ, {y, -ymax, ymax}, 
                    DisplayFunction -> Identity][[1, 1, 1]] /. 
                Line[{a___, {y1_?NonPositive, x1_}, {y2_?NonNegative, x2_}, 
                      b___}] :> {Line[{a, {y1, x1}, {0, 0}}], 
                    Line[{{0, 0}, {y2, x2}, 
                        b}]}, {3}]; 
u = 
            Min[xmax, ymax]; 
coeffList = 
            Table[(Sin[r]*u - Cos[r]*u*Log[Cos[r]*u]/λ)/(Cos[r]*
                    u), {r, -Pi/2 + Pi/n, Pi/2 - Pi/n, 
                Pi/n}]; 
orizz = {style, 
                Map[N[sim] . # &, {Line[{{-xmax, 0}, {0, 0}}], 
                    Line[{{0, 0}, {xmax, 0}}]}, {3}]} // 
              Flatten; 
basicTrajs = {Map[sol, coeffList], 
              orizz, equil}; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[basicTrajs, {-1, 1}, {-1, 1}, 
                wrongDirectionQ[mat], arwDnsty, opts], basicTrajs]; 
          Show[{If[shwFldQ, fld, Sequence @@ {}], Graphics[allTrajs]}, 
            PlotRange -> {{-1, 1}, {-1, 1}}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];

Nonvanishing multiple of the identity matrix

LinearPhasePlot[mat_?nonNullMultipleOfIdentity, opts___Rule] := 
        Module[{λ, wrtRptQ, shwFldQ, shwArrsQ, arwDnsty, n, fld, sol, 
            coeffList, basicTrajs, allTrajs}, 
λ = 
            mat[[1, 1]]; 
{wrtRptQ, shwFldQ, shwArrsQ, 
              arwDnsty, 
              n} = {WriteReport, ShowField, CurvesWithArrows,  
                  ArrowDensity, NumberOfTrajectories} /. {opts} /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print[If[λ > 
                    0, "Repulsive", "Attractive"] <> " node with \
coincident eigenvalues"]]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, - .9,  .9}, {y, - .9,  .9}, 
              DisplayFunction -> Identity]; 
sol[
              c_] := {Line[2 {-{Cos[c], Sin[c]}, {0, 0}}], 
              Line[2 {{0, 0}, {Cos[c], 
                      Sin[c]}}]}; 
coeffList = 
            N[Pi*Range[
                    2  n]/(2  n)]; 
basicTrajs = {Map[
                sol, coeffList], equil}; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[basicTrajs, {-1, 1}, {-1, 1}, 
                wrongDirectionQ[mat], arwDnsty, opts], basicTrajs]; 
          Show[{If[shwFldQ, fld, Sequence @@ {}], Graphics[allTrajs]}, 
            PlotRange -> {{-1, 1}, {-1, 1}}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];

Two null eigenvalues but the matrix is not null

LinearPhasePlot[mat_?nonNullNilpotentQ, opts___Rule] := 
        Module[{sim, canon, xmax, ymax, wrtRptQ, shwFldQ, shwArrsQ, arwDnsty, 
            n, style, fld, sol, orizz, provisional, coeffList, basicTrajs, 
            allTrajs}, 
{sim, canon} = 
            JordanDecomposition[mat]; 
{xmax, ymax} = 
            N[inverseRange[sim]]; 
{wrtRptQ, shwFldQ, 
              shwArrsQ, arwDnsty, n, 
              style} = {WriteReport, ShowField, CurvesWithArrows,  
                  ArrowDensity, NumberOfTrajectories, 
                  NotableLineStyle} /. {opts} /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print["Two null eigenvalues"], Null]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, - .9,  .9}, {y, - .9,  .9}, 
              DisplayFunction -> Identity]; 
sol[c_] := 
            Map[N[sim] . # &, {Line[{{-xmax, c}, {0, c}}], 
                Line[{{0, c}, {xmax, 
                      c}}]}, {3}]; 
orizz = {{GrayLevel[1], 
                Thickness[ .04], 
                provisional[(N[sim] . # &) /@ {{-xmax, 0}, {xmax, 
                        0}}]}, {RGBColor[1, 0, 0], Dashing[{ .02,  .02}], 
                  Thickness[ .01], style, 
                  provisional[(N[sim] . # &) /@ {{-xmax, 0}, {xmax, 
                          0}}]} // Flatten}; 
coeffList = 
            ymax*Complement[Range[-n, n], {0}]/
                n; 
basicTrajs = {Map[sol, coeffList], 
              orizz, equil}; 
allTrajs = 
            If[shwArrsQ, 
                clipAndPlaceArrows[basicTrajs, {-1, 1}, {-1, 1}, 
                  wrongDirectionQ[mat], arwDnsty, opts], basicTrajs] /. 
              provisional -> Line; 
          Show[{If[shwFldQ, fld, Sequence @@ {}], Graphics[allTrajs]}, 
            PlotRange -> {{-1, 1}, {-1, 1}}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];

Null matrix

LinearPhasePlot[mat_?nullMatrixQ, 
          opts___Rule] := 
Show[Graphics[equil], 
          PlotRange -> {{-1, 1}, {-1, 1}}, 
          FilterOptions[Plot, 
            opts], 
If[WriteReport /. {opts} /. 
              Options[LinearPhasePlot], Print["Null matrix"], 
            Null]; 
FilterOptions[Plot, 
            Options[LinearPhasePlot]], 
          DisplayFunction -> $DisplayFunction];

Center (purely imaginary conjugated eigenvalues)

LinearPhasePlot[mat_?centerQ, opts___Rule] := 
        Module[{sim, canon, xmax, ymax, wrtRptQ, shwFldQ, shwArrsQ, arwDnsty, 
            n, style, fld, symAxes, θ0, unitCirc, sol, coeffList, 
            basicTrajs, allTrajs, ax}, 
{sim, canon} = 
            realCanonicalForm[mat]; {xmax, ymax} = 
            N[inverseRange[sim]]; 
{wrtRptQ, shwFldQ, 
              shwArrsQ, arwDnsty, n, 
              style} = {WriteReport, ShowField, CurvesWithArrows,  
                  ArrowDensity, 1 + NumberOfTrajectories, 
                  NotableLineStyle} /. {opts} /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print["Center (periodic orbits)"], 
            Null]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, - .9,  .9}, {y, - .9,  .9}, 
              DisplayFunction -> 
                Identity]; 
symAxes = Eigensystem[
                Transpose[sim] . 
                  sim][[2]]; 
θ0 = 
            N[ArcTan @@ symAxes[[1]]]; 
          unitCirc = 
            Line[Table[{Cos[θ], 
                    Sin[θ]}, {θ, θ0, θ0 + 2  Pi - 
                      2  Pi/80. , 2  Pi/80. }]~
                Join~{{Cos[θ0], 
                    Sin[θ0]}}]; 
sol[c_] := 
            Map[c*N[sim] . # &, unitCirc, {2}]; 
          coeffList = 
            Max[xmax, ymax]*Range[n]^2/n^2; 
basicTrajs = 
            Map[sol, coeffList]; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[basicTrajs, {-1, 1}, {-1, 1}, 
                wrongDirectionQ[mat], arwDnsty, opts], basicTrajs]; 
          ax[i_] := {Dashing[{ .02,  .02}], style, 
                Map[N[sim] . # &, 
                  Line[Sqrt[2] 
                      Max[xmax, ymax]/
                        Sqrt[symAxes[[i]] . 
                            symAxes[[i]]] {-symAxes[[i]], 
                        symAxes[[i]]}], {2}]} // 
              Flatten; 
Show[{If[shwFldQ, fld, 
                Sequence @@ {}], Graphics[{allTrajs, ax[1], ax[2], equil}]}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];

Spiral (complex conjugate eigenvalues with nonvanishing real part)

LinearPhasePlot[mat_?spiralQ, opts___Rule] := 
        Module[{sim, canon, xmax, ymax, α, β, wrtRptQ, shwFldQ, 
            shwArrsQ, arwDnsty, n, style, fld, basicSpir, sol, coeffList, 
            basicTrajs, allTrajs, symAxes, 
            ax}, 
{sim, canon} = 
            realCanonicalForm[mat]; {xmax, ymax} = 
            N[inverseRange[
                sim]]; 
{α, β} = {canon[[1, 
                  1]], canon[[2, 1]]}; 
{wrtRptQ, 
              shwFldQ, shwArrsQ, arwDnsty, n, 
              style} = {WriteReport, ShowField, CurvesWithArrows,  
                  ArrowDensity, 2  NumberOfTrajectories, 
                  NotableLineStyle} /. {opts} /. 
              Options[LinearPhasePlot]; 
If[wrtRptQ, 
            Print[If[α > 
                    0, "Repulsive", "Attractive"] <> " spirals"], 
            Null]; 
fld := 
            PlotVectorField[
              mat . {x, y}, {x, - .9,  .9}, {y, - .9,  .9}, 
              DisplayFunction -> 
                Identity]; 
basicSpir = ParametricPlot[
                Exp[α*θ] {Cos[β*θ], 
                    Sin[β*θ]}, 
                Evaluate[{θ, Log[ .01  Min[xmax, ymax]]/α, 
                    Log[Sqrt[xmax^2 + ymax^2]]/α}], 
                DisplayFunction -> Identity][[1, 1, 
                1]]; 
sol[c_] := 
            Map[N[sim] . {{Cos[c], Sin[c]}, {-Sin[c], Cos[c]}} . # &, 
              basicSpir, {2}]; 
          coeffList = Range[n]*2  Pi/n; 
basicTrajs = 
            Map[sol, coeffList]; 
allTrajs = 
            If[shwArrsQ, 
              clipAndPlaceArrows[basicTrajs, {-1, 1}, {-1, 1}, 
                wrongDirectionQ[mat], arwDnsty, opts], basicTrajs]; 
          symAxes = Eigensystem[Transpose[sim] . sim][[2]]; 
          ax[i_] := {Dashing[{ .02,  .02}], style, 
                Map[N[sim] . # &, 
                  Line[Sqrt[2] 
                      Max[xmax, ymax]/
                        Sqrt[symAxes[[i]] . 
                            symAxes[[i]]] {-symAxes[[i]], 
                        symAxes[[i]]}], {2}]} // 
              Flatten; 
Show[{If[shwFldQ, fld, 
                Sequence @@ {}], Graphics[{allTrajs, ax[1], ax[2], equil}]}, 
            FilterOptions[Plot, opts], 
            FilterOptions[Plot, Options[LinearPhasePlot]], 
            DisplayFunction -> $DisplayFunction]];

Poincaré map and fixed point thereof

Default options and error messages

Options[PoincareMap] = {GridPoints -> 5, 
          GridParameterization -> ({##} &)};

Options[PoincareMapPlot] = {GridStyle -> Automatic, 
          OnlyEndGrid -> False, Superimposed -> False, 
          StackVertical -> False, 
          SamePlotRange -> True};

PoincareMap::"BadGridPoints" = "Bad GridPoints"; 
    PoincareMap::"BadParameterization" = "Bad \
GridParameterization";

This utility processes the grid data given by the user into a form more readily usable by the other functions.

\(\(baseGrid[{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
            GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
            GridVMin < gridVMax, 
          opts___?OptionQ] := 
Module[{hm, grdPrmtrzt, 
            ptsU, ptsV, grdPtsU, grdPtsV, 
            r, θ}, 
{hm, 
              grdPrmtrzt} = \({GridPoints, GridParameterization} /. {opts}\) /. 
              Options[PoincareMap]; 
{ptsU, ptsV} = 
            Which[MatchQ[hm, {_Integer?Positive, __Integer?Positive}], hm, 
              MatchQ[hm, _Integer?Positive], {hm, hm}, True, 
              Message[PoincareMap::"BadGridPoints"]; 
              Abort[]]; 
grdPtsU = 
            GridUMin + \((gridUMax - GridUMin)\) 
                Range[0, 
                    ptsU - 1]/\((ptsU - 1)\); 
grdPtsV = 
            GridVMin + \((gridVMax - GridVMin)\) 
                Range[0, 
                    ptsV - 1]/\((ptsV - 
                      1)\); 
Which[
grdPrmtrzt === Polar, 
            grdPrmtrzt = 
              Function[{r, θ}, r {Cos[θ], Sin[θ]}], 
            grdPrmtrzt === Identity, grdPrmtrzt = \(({##} &)\), 
            Not[And[Map[NumericQ, grdPrmtrzt[GridUMin, GridVMin]]]], 
            Message[PoincareMapPlot::"BadParameterization"; 
              Abort[]]]; 
{grdPtsU, grdPtsV, 
            grdPrmtrzt}];\)\)

This utility plots a parametric grid.

gridPlot[grdPtsU_List, grdPtsV_List, grdPrmtrzt_, opts___?OptionQ] := 
        Module[{grdstl, ptsU, ptsV, stlU, stlV, a, b, trajU, 
            trajV}, 
grdstl = GridStyle /. {opts} /. 
              Options[PoincareMapPlot]; 
{ptsU, ptsV} = 
            Map[Length, {grdPtsU, grdPtsV}]; 
stlU = 
            Which[grdstl === Automatic, 
              Map[{Hue[#]} &, Range[0, ptsV]/(2  ptsV)], 
              grdstl == None, Array[{} &, grdPtsV], True, 
              grdstl[[1]]]; 
          stlV = Which[grdstl === Automatic, 
              Map[{Hue[#]} &, -Range[1, ptsU]/(2  ptsU)], 
              grdstl == None, Array[{} &, grdPtsU], True, 
              grdstl[[2]]]; 
trajU = 
            ParametricPlot[
              Evaluate[Map[grdPrmtrzt[a, #] &, grdPtsV]], {a, Min[grdPtsU], 
                Max[grdPtsU]}, DisplayFunction -> Identity, 
              PlotStyle -> stlU, 
              Evaluate[
                FilterOptions[ParametricPlot, 
                  opts]]]; 
trajV = 
            ParametricPlot[
              Evaluate[Map[grdPrmtrzt[#, b] &, grdPtsU]], {b, Min[grdPtsV], 
                Max[grdPtsV]}, DisplayFunction -> Identity, 
              PlotStyle -> stlV, 
              Evaluate[
                FilterOptions[ParametricPlot, 
                  opts]]]; 
Show[{trajU, trajV}, 
            DisplayFunction -> Identity, FilterOptions[Graphics, opts], 
            FilterOptions[Graphics, 
              Sequence @@ Options[PoincareMapPlot]]]];

This function calculates the Poincaré map on a grid of points and interpolates between the values.

PoincareMap[data___] = internalPoincareMap[data];

internalPoincareMap[{{x_Symbol, y_Symbol}, {xprime_, 
              yprime_}}, {t_Symbol, t0_?NumericQ, 
            tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
              GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
              GridVMin < gridVMax}, 
          opts___?OptionQ] := internalPoincareMap[{{x, y}, {xprime, 
                yprime}}, {t, t0, 
              tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, opts] = 
          Module[{grdPtsU, grdPtsV, grdPrmtrzt, strtPtsGrid, endPt, x0, y0, 
              u, v, endPtsGrid, 
              pairOfInterp}, 
{grdPtsU, grdPtsV, 
                grdPrmtrzt} = 
              baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
                opts]; 
strtPtsGrid = 
              Apply[grdPrmtrzt, 
                Outer[List, grdPtsU, grdPtsV], {2}]; 
endPt[
                x0_, y0_] := {u[tMax], 
                  v[tMax]} /. NDSolve[{u'[t] == 
                          xprime, v'[t] == yprime, u[0] == x0, 
                        v[0] == y0} /. {x -> u[t], 
                        y -> v[t]}, {u, v}, {t, t0, tMax}, 
                    FilterOptions[NDSolve, 
                      opts]][[1]]; 
endPtsGrid = 
              Apply[endPt, 
                strtPtsGrid, {2}]; 
pairOfInterp = 
              Map[ListInterpolation[#, {{GridUMin, gridUMax}, {GridVMin, 
                        gridVMax}}] &, 
                Transpose[
                  endPtsGrid, {2, 3, 1}]]; 
Function[
              Evaluate[Through[pairOfInterp[#1, #2]]]]];

This function plots the Poincaré map.

PoincareMapPlot[{{x_Symbol, y_Symbol}, {xprime_, 
              yprime_}}, {t_Symbol, t0_?NumericQ, 
            tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
              GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
              GridVMin < gridVMax}, 
          opts___?OptionQ] := 
Module[{onlyEndQ, superImpQ, 
            stkVertQ, samePR, grdPtsU, grdPtsV, grdPrmtrzt, baseGridPlot, 
            poinc, poincPlot, 
            plRng}, 
{onlyEndQ, superImpQ, stkVertQ, 
              samePR} = {OnlyEndGrid, Superimposed, StackVertical, 
                  SamePlotRange} /. {opts} /. 
              Options[PoincareMapPlot]; 
{grdPtsU, grdPtsV, 
              grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
baseGridPlot := 
            gridPlot[grdPtsU, grdPtsV, grdPrmtrzt, 
              DisplayFunction -> Identity, 
              opts]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
poincPlot = 
            gridPlot[grdPtsU, grdPtsV, poinc, 
              DisplayFunction -> Identity, 
              opts]; 
plRng := 
            If[samePR == True, 
              If[$VersionNumber < 4, 
                PlotRange -> 
                  FullOptions[Show[{baseGridPlot, poincPlot}], PlotRange], 
                AbsoluteOptions[Show[{baseGridPlot, poincPlot}], PlotRange]], 
              Sequence @@ {}]; 
Show[
Which[
onlyEndQ === 
                True, poincPlot, 
superImpQ, {baseGridPlot, 
                poincPlot}, 
stkVertQ, 
              GraphicsArray[{{Show[baseGridPlot, plRng]}, {Show[poincPlot, 
                      plRng]}}], 
True, 
              GraphicsArray[{Show[baseGridPlot, plRng], 
                  Show[poincPlot, plRng]}]], 
FilterOptions[
              Graphics, opts], 
            DisplayFunction -> $DisplayFunction]];

Let φ(x, y) be the Poincaré map at some time for some system. Consider the function ψ(x, 
          y) := ∥φ(x, y) - (x, 
              y) (∥ ^2). The following function makes a DensityPlot of Min[ψ(x, y), ϵ] for a give tolerance ϵ. This way the point of the plot that are gray to black are candidates for being a fixed point for φ.

PoincareMapFixedPointDensityPlot[{{x_Symbol, y_Symbol}, {xprime_, 
                yprime_}}, {t_Symbol, t0_?NumericQ, 
              tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
                GridUMin < gridUMax, {GridVMin_?NumericQ, 
                  gridVMax_?NumericQ} /; GridVMin < gridVMax}, 
            Optional[tolerance_, (1/100)], opts___?OptionQ] /; 
          tolerance > 0 := 
Module[{grdPtsU, grdPtsV, 
            grdPrmtrzt, poinc, diff, u, 
            v}, 
{grdPtsU, grdPtsV, grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
rectGrid = 
            Last[baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
                GridParameterization -> Identity, 
                opts]]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
diff[u_, v_] := (# . # &)[
              poinc[u, v] - 
                grdPrmtrzt[u, v]]; 
Show[{DensityPlot[
                Min[diff[u, v], tolerance], {u, GridUMin, gridUMax}, {v, 
                  GridVMin, gridVMax}, DisplayFunction -> Identity, 
                Evaluate[FilterOptions[DensityPlot, opts]], 
                Mesh -> False], 
              gridPlot[grdPtsU, grdPtsV, rectGrid, 
                DisplayFunction -> Identity, PlotPoints -> 2]}, 
            FilterOptions[Graphics, opts], 
            DisplayFunction -> $DisplayFunction, 
            PlotRange -> All]];

The following function makes a 3D plot of the function ψ  instead of a DensityPlot as above.

PoincareMapFixedPointPlot3D[{{x_Symbol, y_Symbol}, {xprime_, 
              yprime_}}, {t_Symbol, t0_?NumericQ, 
            tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
              GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
              GridVMin < gridVMax}, 
          opts___?OptionQ] := 
Module[{grdPtsU, grdPtsV, 
            grdPrmtrzt, poinc, diff, u, 
            v}, 
{grdPtsU, grdPtsV, grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
diff[u_, v_] := (# . # &)[
              poinc[u, v] - grdPrmtrzt[u, v]]; 
          ParametricPlot3D[
            Join[grdPrmtrzt[u, v], {diff[u, v]}], {u, GridUMin, 
              gridUMax}, {v, GridVMin, gridVMax}, 
            Evaluate[FilterOptions[ParametricPlot3D, opts]]]];

The following function tries to find a fixed point of the Poincaré map by minimizing the function ψ defined above near a point (u0,v0) given by the user.

FindPoincareMapFixedPoint[{{x_Symbol, y_Symbol}, {xprime_, 
                yprime_}}, {t_Symbol, t0_?NumericQ, 
              tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
                GridUMin < gridUMax, {GridVMin_?NumericQ, 
                  gridVMax_?NumericQ} /; GridVMin < gridVMax}, {u_Symbol, 
              u0__?NumericQ}, {v_Symbol, v0__?NumericQ}, 
            opts___?OptionQ] /; (GridUMin ≤ 
                First[{u0}] ≤ Last[{u0}] ≤ gridUMax && 
              GridVMin ≤ First[{v0}] ≤ 
                Last[{v0}] ≤ gridVMax) := 
        Module[{grdPtsU, grdPtsV, grdPrmtrzt, poinc, diff, a, 
            b}, 
{grdPtsU, grdPtsV, grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
diff[a_, b_] := (# . # &)[
              poinc[a, b] - grdPrmtrzt[a, b]]; 
FindMinimum[
            diff[u, v], {u, u0}, {v, v0}, 
            Evaluate[FilterOptions[FindMinimum, opts]]]];

Tail end

End[];

Protect["`*"];

EndPackage[];

On[General::"spell1"]

On[General::"spell"]

Tutorial examples

Evaluate the input cells below to see the package in action and learn its usage.

OrientedParametricPlot

OrientedParametricPlot works just like ParametricPlot, but places arrow heads along the line:

OrientedParametricPlot[{Cos[x], Sin[x]}, {x, 0, 3  Pi/2}, 
        HowManyArrows -> 3, AspectRatio -> Automatic];

You can dictate the exact placement of the arrows along the curve with ArrowPositions. Values must be between 0 and 1, and will place arrow heads proportionally with respect to arc length (allow for some distoriton when AspectRatio is not set to Automatic).

OrientedParametricPlot[{Cos[x], Sin[x]}, {x, 0, 3  Pi/2}, 
        ArrowPositions -> {1/4, 1/2, 1}, 
        AspectRatio -> Automatic];

A more complicated example:

OrientedParametricPlot[{Cos[x] + Sin[3  x]/3, Sin[x]}, {x, 0, 
          3  Pi/2}, HowManyArrows -> 3, 
        AspectRatio -> Automatic];

You can reverse the direction of the arrows without changing parameterization:

OrientedParametricPlot[{Cos[x], Sin[x]}, {x, 0, 3  Pi/2}, 
        HowManyArrows -> 2, 
        ReverseDirectionQ -> True];

More of one curve on the same plot is allowed, just as in ParametricPlot:

OrientedParametricPlot[{{Cos[x], Sin[x]}, {2  Cos[x], 
            2  Sin[x]}}, {x, 0, 3  Pi/2}, HowManyArrows -> 2, 
        PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}}, 
        AspectRatio -> Automatic, 
        Prolog -> {RGBColor[0, 0, 1], PointSize[ .04], 
            Point[{0, 0}]}];

Trig Interpolation curves

I have written some utilities to help me draw pictures of "generic-looking" smooth functions or curves in the plane. I was not satisfied with polynomial interpolation or splines, because that way you can just spot the junctions, even if they are technically smooth. So I turned to trigonometric functions instead of polynomials to see what happened.

One-variable interpolation

TrigInterpolation is an alternative to Interpolation: it tries to find a finite-Fourier-series functions that fits your data points. With TrigInterpolation you get extrapolation as well, because the function you get is defined everywhere. The basic syntax is the same as for Interpolation:

TrigInterpolation[{{0, 1}, {1, 1}, {2, 0}}, x]

TrigInterpolation[{{1, 1, 0}, {2, 0}, {3, 4}}, t]

When you have more than a few conditions, consider using floating-point input, so that the computation is faster and you avoid unwieldy formulas. Compare the following two cells:

data = {{1, 1, 0}, {2, 0}, {3, 4, 0}, {5, 
            1}};
TrigInterpolation[data, t]

TrigInterpolation[data // N, t]

If you insist on exact input, with the option EvaluateOnCoefficients you can assign a function to apply to the coefficients of the trigonometric sum. For example you can FullSimplify the coefficients, without affecting the trig functions:

TrigInterpolation[{{1, 1, 0}, {2, 0}, {3, 4, 0}, {5, 1}},t,EvaluateOnCoefficients -> FullSimplify]

The function TrigInterpolationPlot draws for you the plot of the trig interpolation. The syntax is TrigInterpolationPlot[data,xmin,xmax,options]. The output is a list: the first element is the graph and the last element is the output of TrigInterpolation (in pure function format):

TrigInterpolationPlot[{{1, 1}, {2, 0}, {3, 4}} // N, 0, 4]

With the option ShowListPlot you can overimpose the interpolation points to the graph of the intepolating function:

TrigInterpolationPlot[{{1, 1}, {2, 0}, {3, 4}} // N,0,4,ShowListPlot -> True]

You can control the style of the interpolation points with ListPlotStyle:

TrigInterpolationPlot[{{1, 1}, {2, 0}, {3, 4}} // N,0,4,ShowListPlot -> True,ListPlotStyle -> {RGBColor[1, 0, 0], 
              PointSize[ .03]}]

The period of the interpolating trig function is by default 3 times the range of the independent variable for the data points. You can choose your own period with Period.

TrigInterpolationPlot[{{1, 1}, {2, 0}, {3, 4}} // N,0,4,ShowListPlot -> True,Period -> 2.5]

TrigInterpolationPlot[{{1, 1}, {2, 0}, {3, 4}} // N,0,4,ShowListPlot -> True,Period -> 10]

Beware that some period values may be incompatible with the data.

TrigInterpolationPlot[{{1, 1}, {2, 2}}, 0, 4, 
      Period -> 1]

Curve interpolation

OpenCurveInterpolationPlot[{{{0, 0}, {0, 0}, {0, 1}}, {1, 1}, {0, 1}} // 
        N, ShowListPlot -> True, Axes -> False, 
      ShowTimesList -> True]

PolynomialLevelSetInterpolation[{{0, 0}, {1, 0}, {1, 1}, {0, 2}}, t, x, 
      y, ShowTimesList -> True]

ClosedCurveInterpolation[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}, 
      t]

ClosedCurveInterpolation[{{{0, 0}, {0, 0}, {0, 1}}, {1, 1}, {0, 1}} // 
        N, t, ShowTimesList -> True]

ClosedCurveInterpolation[{{{0, 0}, {1, 0}}, {1, 0}, {1, 1}, {0, 1}}, t, 
      EvaluateOnCoefficients -> RootReduce]

ClosedCurveInterpolationPlot[{{{0, 0}, {0, 0}, {0, 1}}, {1, 1}, {0, 
            1}} // N, ShowListPlot -> True, 
      Axes -> False]

ClosedCurveInterpolationPlot[{{{0, 0}, {1/2, -1/2}}, {1, 1}, {0, 
            1}} // N, Axes -> True, ShowListPlot -> True, 
      ListPlotStyle -> PointSize[ .03]]

ClosedCurveInterpolationPlot[{{{0, 0}}, {1, 1}, {0, 1}} // N, 
      Axes -> True, ShowListPlot -> True]

OrientedClosedCurveInterpolationPlot[{{{0, 0}, {0, 0}, {0, 1}}, {1, 
            1}, {0, 1}} // N, ShowListPlot -> True, Axes -> False, 
      ShowTimesList -> True]

FilledClosedCurveInterpolationPlot[{{{0, 0}, {0, 0}, {0, 1}}, {1, 
            1}, {0, 1}} // N, ShowListPlot -> True, 
      PlotStyle -> Dashing[{ .01,  .02}], Fills -> Red]

Level set interpolation

ClosedCurveInterpolationPlot[{{0, 0}, {2, 0}, {1, 1}, {0, 
            1}}];

PolynomialLevelSetInterpolation[{{0, 0}, {2, 0}, {1, 1}, {0, 1}}, t, x, 
      y]

ContourPlot[
        ArcTan[2 x - (x^2) + 4 y - x y - 8 (y^2) + 5 (y^3) - 
            (y^4)], {x, -1, 3}, {y, -1, 2}, 
        Contours -> {-1/2, 0, 1/2}];

NDSolvePlot

A first-order planar system:

NDSolvePlot[{x'[t] == x[t]/10 - y[t], y'[t] == 
          y[t]/10 + x[t], x[0] == 1, y[0] == 0}, {x, y}, {t, 0, 
        10}, HowManyArrows -> 3]

A second-order planar system:

NDSolvePlot[{x''[t] == -x[t] (1 + x[t]^2), y''[
            t] == -y[t] (1 + x[t]^2 - x[t]^4), x[0] == 1, 
        y[0] == 0, x'[0] == 0, y'[0] == 1}, {x, 
        y}, {t, 0, 10}, HowManyArrows -> 6, 
      AspectRatio -> Automatic]

The output is a list {NDSolve output, Graphics output}.

A second-order scalar equation (the default for scalar equations is without arrows):

NDSolvePlot[{x''[t] + x'[t]/2 == -x[t], 
        x[0] == 1, x'[0] == 0}, x, {t, 0, 10}, 
      PlotRange -> All]

Same as above, but with arrows:

NDSolvePlot[{x''[t] + x'[t]/2 == -x[t], 
          x[0] == 1, x'[0] == 0}, x, {t, 0, 10}, 
        PlotRange -> All, HowManyArrows -> 3];

You can give the same options as with OrientedParametricPlot:

NDSolvePlot[{x'[t] == y[t], y'[t] == 1, 
          x[0] == 0, y[0] == 0}, {x, y}, {t, 0, 1}, 
        HowManyArrows -> 2, PlotStyle -> {RGBColor[1, 0, 0]}, 
        Prolog -> {RGBColor[1, 0, 0], PointSize[ .03], 
            Point[{0, 0}]}];

HamiltonianPlot

HamiltonianPlot draws a phase portrait of a Hamiltonian system, based on ContourPlot with arrow heads added:

HamiltonianPlot[
        x^2 + y^2, {x, -1, 1}, {y, -1, 1}];

You can add the direction field:

HamiltonianPlot[x^2 + y^2, {x, -1, 1}, {y, -1, 1}, 
        ShowField -> True];

The pendulum system:

HamiltonianPlot[-Cos[x] + y^2/2, {x, -Pi, Pi}, {y, -Pi, 
          Pi}];

Increasing the option PlotPoints you can make the curves less jagged:

HamiltonianPlot[-Cos[x] + y^2/2, {x, -Pi, Pi}, {y, -Pi, 
          Pi}, PlotPoints -> 30];

PhasePlot

PhasePlot draws trajectories of a planar auronomous system of ordinary differential equations, placing arrows along the curves:

PhasePlot[{1/(1 + x^2), -y + x}, {x, -1, 1}, {y, -1, 
          1}, {0, 2}];

LinearPhasePlot

LinearPhasePlot takes a two-by-two real matrix as argument, plus options.The matrix must be made of exact numbers, not of floating-point numbers.

The function LinearPhasePlot accepts all Plot options, but it ignores PlotRange.

With the option NotableLineStyle you can control how the following special lines are drawn: the eigenspaces when the eigenvalues are real, and the symmetry axes of the ellipse image of the unit circle when the eigenvalues are complex.

With the option NumberOfTrajectories you can control how many trajectories are drawn. The default value is 5, which was found to be reasonable. If you want less trajectories try a value of 4 or 3; if you want more try 6 or 7.

With the option WriteReport set to True (the default) you get a classification of the kind of equilibrium the origin is.

Try out the following examples:

Repulsive node:

LinearPhasePlot[{{1, 1/2}, {1/5, 2}}];

Same as above, with the field:

LinearPhasePlot[{{1, 1/2}, {1/5, 2}}, 
        ShowField -> True];

Same as above, with eigenspaces colored greeen instead of red:

LinearPhasePlot[{{1, 1/2}, {1/5, 2}}, 
        NotableLineStyle -> RGBColor[0, 1, 0]];

By turning off the insertion of arrow along the curves (as distinct from the grid-arrows of the field) the computation is faster:

LinearPhasePlot[{{1, 1/2}, {1/5, 2}}, 
        CurvesWithArrows -> False];

Saddle:

LinearPhasePlot[{{-1, 0}, {1, 1}}];

Saddle with the field:

LinearPhasePlot[{{-1, 0}, {1, 1}}, 
        ShowField -> True];

One null eigenvalue, the other one positive:

LinearPhasePlot[{{1, 2}, {3, 4}} . {{1, 0}, {0, 0}} . 
          Inverse[{{1, 2}, {3, 4}}]];

Improper node:

LinearPhasePlot[{{2, 1}, {-1, 4}}];

Identity matrix:

LinearPhasePlot[-{{1, 0}, {0, 1}}];

Two null eigenvalues:

LinearPhasePlot[{{1, 2}, {3, 4}} . {{0, 1}, {0, 0}} . 
          Inverse[{{1, 2}, {3, 4}}]];

Null matrix:

LinearPhasePlot[{{0, 0}, {0, 0}}];

Center:

LinearPhasePlot[{{1, 2}, {3, 1}} . {{0, 1}, {-1, 0}} . 
          Inverse[{{1, 2}, {3, 1}}]];

Spiral (it seems to take more time to compute than the other examples):

LinearPhasePlot[{{1/2, -1}, {1, 1/2}}];

Poincaré maps and their fixed points

Simple examples of Poincaré map plots

PoincareMapPlot takes a system of two differential equations of the form {({{x' = a(t, x, y)}, {y' = b(t, x, y)}})and draws the initial and final position of a user-specified grid under the Poincaré map at a given time. The system is specified with the syntax {{x, y}, {a, b}}, where x, y are variable names and a, b the right-hand sides:

field = {{x, y}, {1/(1 + x^2), -y + x}};

The grid is given as {{xmin, xmax}, {ymin, ymax}}.

grid = {{0, 1}, {0, 1}};

Here is the graph of the Poincaré map at time 1/2 with all default settings. On the left is the original grid at time 0, and on the right where the grid has drifted at time 1/2

PoincareMapPlot[field, {t, 0, 1/2}, grid];

With the option GridPoints you can control the number of lines in the grid. The default is 5 in each direction:

PoincareMapPlot[field, {t, 0, 1/2}, grid, 
        GridPoints -> {3, 5}];

The grid does not need to be cartesian. You can assign a polar grid, for example:

PoincareMapPlot[field, {t, 0, 1/2}, {{0, 1}, {0, 2  Pi}}, 
        GridParameterization -> Polar, AspectRatio -> Automatic, 
        GridPoints -> {5, 10}];

You can have an arbitrarily parameterized grid:

PoincareMapPlot[field, {t, 0, 1/2}, {{0, 1}, {0, 1}}, 
        GridParameterization -> Function[{a, b}, {a + b, a - b}], 
        AspectRatio -> Automatic, 
        GridPoints -> {5, 10}];

You can stack the two plots vertically:

PoincareMapPlot[field, {t, 0, 1/2}, {{0, 1}, {0, 1}}, 
        GridParameterization -> Function[{a, b}, {a + b, a - b}], 
        AspectRatio -> Automatic, GridPoints -> {5, 10}, 
        StackVertical -> True];

You can also decide to draw the final grid only:

PoincareMapPlot[field, {t, 0, 1/2}, {{0, 1}, {0, 1}}, 
        GridParameterization -> Function[{a, b}, {a + b, a - b}], 
        AspectRatio -> Automatic, GridPoints -> {5, 10}, 
        OnlyEndGrid -> True, AxesOrigin -> {0, 0}];

You can also have the initial and final grid superimposed:

PoincareMapPlot[field, {t, 0, 1/2}, {{0, 1}, {0, 1}}, 
        GridParameterization -> Function[{a, b}, {a + b, a - b}], 
        AspectRatio -> Automatic, GridPoints -> {5, 10}, 
        Superimposed -> True];

By default the initial and final grid are drawn with the same PlotRange. You can override this behaviour and let each plot go on its own with the option SamePlotRange->False.

PoincareMapPlot[field, {t, 0, 1}, {{0, 1}, {0, 1}}, 
        GridParameterization -> Function[{a, b}, {a + b, a - b}], 
        GridPoints -> {5, 10}, 
        SamePlotRange -> False];

Poincaré map of x'' (t) = -f(x(t)), and search for periodic solutions of a given period.

Take the following function of one variable

f = x/(1 + x^2) + 3  x (x^2/(1 + x^2));

Here in red is the graph of f. Near x = 0 f(x) ≈ x, while f(x) ≈ 3  x for large x.

Plot[{x, f, 3  x}, {x, 0, 4}, 
        PlotStyle -> {{}, {Thickness[ .01], RGBColor[1, 0, 0]}, {}}, 
        AspectRatio -> 1];

Solutions to x'' (t) = -f(x(t)) that start with small values of (x^2) + (((x'))^2) will be periodic of period approximately  2  π ≈ 6.28, and those with (x^2) + (((x'))^2)large will have period approximately  ((2 π)/(Sqrt[3])) ≈ 3.63. We will try to locate a solution with an intermediate period, namely 5.

The program wants the second-order differential equation x'' (t) = -f(x(t)) to be framed as a system of first-order equation of the general form {({{x' = a(t, x, y)}, {y' = b(t, x, y)}}), and the data to be arranged as {{x, y}, {a(t, x, y), b(t, x, y)}}:

field = {{x, y}, {y, -f}};

Let us take a polar grid with 0 ≤ r ≤ 2, \(TraditionalForm\`\(-π\) ≤ θ ≤ π\):

grid = {{0, 2}, {-Pi, Pi}};

Here is the graph of the Poincaré map at time 5:

PoincareMapPlot[field, {t, 0, 5}, grid, GridPoints -> {10, 41}, 
        GridParameterization -> Polar, 
        AspectRatio -> Automatic];

The plot on the left is the starting grid, while on the right you see the position of the grid lines after the map is applied. If for any reason you only want the plot on the right just add the option OnlyEndGrid->True:

PoincareMapPlot[field, {t, 0, 5}, grid, GridPoints -> {10, 41}, 
        GridParameterization -> Polar, AspectRatio -> Automatic, 
        OnlyEndGrid -> True];

You can also have the two grids on top of each other with the option Superimposed->True. You had better do that with less dense grids than the present one. You will se an example later.

If  ϕ(x, y) the Poincaré map at time 5, consider the real function d(x, 
          y) := ∥ϕ(x, y) - (x, 
              y) (∥ ^2). A very small value for d means a solution pretty much periodic with period 5. Here is a 3D plot of  d(x, y):

PoincareMapFixedPointPlot3D[field, {t, 0, 5}, grid, 
        GridPoints -> {10, 41}, PlotPoints -> {20, 41}, 
        GridParameterization -> Polar, 
        AxesLabel -> {"x", "y", "z"}, 
        PlotRange -> {0, 1}];

With a density plot of Min[d(x, y), ϵ] the regions where d(x, y)is smaller than ϵ will stand out as grey or black against a white background. In the following picture ϵ = 1/100. The parameters on the two axes are not x, y but r, θ:

PoincareMapFixedPointDensityPlot[field, {t, 0, 5}, grid, 1/100, 
        GridPoints -> {20, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

You will see two distinct grey-black stripes: one around  r = 0, which corresponds to the trivial periodic solution x(t) ≡ y(t) ≡ 0, whilst the other seems to be the graph of a function from θ to r. Let us zoom around the gray point r =  .75, θ = 0:

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 5}, {{ .65,  .85}, {- .2,  .2}}, 1/1000, 
        GridPoints -> {20, 21}, PlotPoints -> {20, 20}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

Let us close in on the fixed point by minimizing the function  d near  r =  .756, θ = 0 with the built-in minimization routines:

{squareDistance, initialData} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 
          5}, {{ .65,  .85}, {- .1,  .1}}, {r,  .76}, {θ, 0}, 
        GridPoints -> {40, 11}, 
        GridParameterization -> Polar]

Here is a plot of the solution corresponding the fixed point found above, seen as a function t ↦ x(t):

{{solutionX}, graphOfX} = 
        NDSolvePlot[{x''[t] == (-f /. x -> x[t]), 
            x[0] == r*Cos[θ] /. initialData, x'[0] == 
              0}, x, {t, 0, 5}];

A check to see how much the values at t = 0 and t = 5 are close to each other:

{x[5] - x[0], x'[5] - x'[0]} /. solutionX

A graph of our periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

{{solutionXY}, phaseSpacePlot} = 
        NDSolvePlot[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == r*Cos[θ] /. initialData, 
            y[0] == 0}, {x, y}, {t, 0, 5}, 
          AspectRatio -> Automatic, HowManyArrows -> 3, 
          PlotStyle -> {Thickness[ .007], 
              Dashing[{ .02,  .04}]}];

The following picture superimposes the 5-periodic solution to the Poincare map plot:

Show[{PoincareMapPlot[field, {t, 0, 5}, grid, 
            GridPoints -> {10, 41}, GridParameterization -> Polar, 
            OnlyEndGrid -> True, AspectRatio -> Automatic, 
            DisplayFunction -> Identity], phaseSpacePlot}, 
        DisplayFunction -> $DisplayFunction];

A zoom of Poincaré map around the starting point of the 5-periodic solution, with the initial and the final grids superimposed:

Show[{phaseSpacePlot, 
          PoincareMapPlot[
            field, {t, 0, 
              5}, {0.756650394233193` + {- .1,  .1}, {- .05,  .05}}, 
            GridPoints -> {7, 5}, GridParameterization -> Polar, 
            Superimposed -> True, AspectRatio -> Automatic, 
            DisplayFunction -> Identity]}, 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> {{ .62,  .86}, {- .282,  .22}}, 
        Frame -> True, Axes -> False];

Animated Poincaré maps

Plots of Poincaré maps at a sequence of times can be shown as a movie. When all the frames have been produced, double-click the first cell in the sequence to start the animation

f = x/(1 + x^2) + 3  x (x^2/(1 + x^2));

field = {{x, y}, {y, -f}};

grid = {{1/4, 2}, {-1, 1}};

filmFrame[instant_] := 
        PoincareMapPlot[field, {t, 0, instant}, grid, 
          OnlyEndGrid -> True, GridPoints -> {5, 5}, 
          GridParameterization -> Polar, AspectRatio -> Automatic, 
          PlotRange -> {{-3, 3}, {-3, 3}}, 
          AxesOrigin -> {0, 0}];

Do[filmFrame[τ], {τ, 0, 5, 1/5}]

Poincaré map of x'' (t) = -f(t, x(t)), and search for periodic solutions of period 2  π.

Take a function f(t, x) that is 2  π-periodic with respect to t and that is asymptotic to x/2 as x -> 0 and to 3  x as x -> +∞:

f = 
        x/(2 + x^2) + 3  x (x^2/(1 + x^2)) + 
          x^2*Cos[t]/(3  x^4 + 2);

Solutions to x'' (t) = -f(t, x(t)) that start with small values of (x^2) + (((x'))^2) will be periodic of period approximately  (Sqrt[8]) π ≈ 8.9, and those with (x^2) + (((x'))^2)large will have period approximately  ((2 π)/(Sqrt[3])) ≈ 3.63. We will try to locate solutions which are periodic with the same t-period as f(t, x), that is, 2  π-periodic solutions. Any such solution will correspond to a fixed point for the Poincaré map at time 2  π.

field = {{x, y}, {y, -f}};

grid = {{0, 2}, {-Pi, Pi}};

Here is a plot of Poincaré map:

PoincareMapPlot[field, {t, 0, 2  Pi}, grid, 
        GridPoints -> {20, 41}, GridParameterization -> Polar, 
        AspectRatio -> Automatic];

If  ϕ(x, y) the Poincaré map at time  2  π, consider the real function d(x, 
          y) := ∥ϕ(x, y) - (x, 
              y) (∥ ^2). A very small value for d means a solution pretty much periodic with period 2  π. With a density plot of Min[d(x, y), ϵ] the regions where d(x, y)is smaller than ϵ will stand out as grey or black against a white background. In the following picture ϵ = 1/1000. The parameters on the two axes are not x, y but r, θ:

PoincareMapFixedPointDensityPlot[field, {t, 0, 2  Pi}, grid, 1/1000, 
        GridPoints -> {40, 81}, PlotPoints -> {100, 100}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

We make out one grey stripe plus two grey spots. The stripe around r = 0 corresponds to the trivial null solution. The two grey spots are approximately at the points ( .45, 0) and ( .7, ±π).

A zoom around (.45,0):

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 2  Pi}, {{ .42,  .48}, {- .5,  .5}}, 1/10000, 
        GridPoints -> {40, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

Minimizing d(x, y) with numerical methods we close in on the fixed point:

{squareDistance1, initialData1} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 
          2  Pi}, {{ .42,  .48}, {- .5,  .5}}, {r,  .445,  .446}, {θ, - .01,  .01}, GridPoints -> {40, 11}, 
        GridParameterization -> Polar]

Here is a plot of the solution corresponding the fixed point found above, seen as a function t ↦ x(t):

{{solutionX1}, graphOfX1} = 
        NDSolvePlot[{x''[t] == (-f /. x -> x[t]), 
            x[0] == r*Cos[θ] /. initialData1, x'[0] == 
              0}, x, {t, 0, 2  Pi}, 
          Ticks -> {{Pi/2, Pi, 3  Pi/2, 2  Pi}, 
              Automatic}];

A check to see how much the values at t = 0 and t = 2  π are close to each other:

{x[2  Pi] - x[0], x'[2  Pi] - x'[0]} /. solutionX1

A graph of our periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

{{solutionXY1}, phaseSpacePlot1} = 
        NDSolvePlot[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == r*Cos[θ] /. initialData1, 
            y[0] == 0}, {x, y}, {t, 0, 2  Pi}, 
          AspectRatio -> Automatic, 
          HowManyArrows -> 3];

A zoom of Poincaré map around the starting point of this 2  π-periodic solution, with the initial and the final grids superimposed:

PoincareMapPlot[
        field, {t, 0, 
          2  Pi}, {{0.44404 -  .1, 0.44404 +  .1}, {- .1,  .1}}, 
        GridPoints -> {5, 5}, GridParameterization -> Polar, 
        AspectRatio -> Automatic, Superimposed -> True, 
        Frame -> True, Axes -> False];

A zoom of the density plot around the second grey spot, the one near ( .7, ±π):

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 2  Pi}, {{ .6,  .8}, {Pi -  .1, Pi +  .1}}, 1/10000, 
        GridPoints -> {40, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

Minimize d(x, y) in that spot:

{squareDistance2, initialData2} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 2  Pi}, {{ .6,  .8}, {3, 3.2}}, {r,  .72}, {θ, 
          Pi}, GridPoints -> {40, 11}, 
        GridParameterization -> Polar]

Here is a plot of the solution corresponding the second fixed point, seen as a function t ↦ x(t):

{{solutionX2}, graphOfX2} = 
        NDSolvePlot[{x''[t] == (-f /. x -> x[t]), 
            x[0] == r*Cos[θ] /. initialData2, x'[0] == 
              0}, x, {t, 0, 2  Pi}];

A check to see how much the values at t = 0 and t = 2  π are close to each other:

{x[5] - x[0], x'[5] - x'[0]} /. solutionX2

A graph of our second periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

{{solutionXY2}, phaseSpacePlot2} = 
        NDSolvePlot[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == r*Cos[θ] /. initialData2, 
            y[0] == 0}, {x, y}, {t, 0, 2  Pi}, 
          AspectRatio -> Automatic, 
          HowManyArrows -> 3];

A zoom of Poincaré map around the starting point of the second 2  π-periodic solution, with the initial and the final grids superimposed:

PoincareMapPlot[
        field, {t, 0, 
          2  Pi}, {{0.70665 -  .1, 0.70665 +  .1}, {Pi -  .1, Pi +  .1}}, 
        GridPoints -> {5, 5}, GridParameterization -> Polar, 
        AspectRatio -> Automatic, Superimposed -> True, 
        Frame -> True, Axes -> False];

Poincaré map of x'' (t) = -f(x(t)) for which there is an invariant annulus

Consider this function of x:

f = 2  x^3;

The solutions of x'' = -f(x) starting with small (x^2) + (((x'))^2) will have large period, while the ones starting with large (x^2) + (((x'))^2) will have small period. We set out to find periodic solutions with period 4.

The orbits of this system, seen as {({{x' = y}, {y' = -f(x)}}) lie on the level sets of the function (x, y) ↦ (x^4) +(y^2), which look somewhat like ellipses. What follows is a pseudo-polar parameterization of the plane where  θ is the usual argument but r = (Sqrt[(x^4) + (y^2)]). The advantage of this parameterization is that the sets r = constant coincide with the level sets of (x, y) ↦ (x^4) +(y^2). The main parameterization has a singularity along the y-axis, so that we need a separate formula for that case.

myParameterization[
          r_?NumericQ, θ_?NumericQ /; 
            Abs[Cos[θ]] < $MachineEpsilon] = {0, 
          r*Sign[Sin[θ]]};
myParameterization[
          r_?NumericQ, θ_?NumericQ /; 
            Abs[Cos[θ]] ≥ $MachineEpsilon] = 
        Sign[Cos[θ]] (Sqrt[-(1/2) Tan[θ]^2 + 
                (1/2) (Sqrt[4 (r^2) + Tan[θ]^4])]) {1, 
            Tan[θ]};

The data structure for the differential system:

field = {{x, y}, {y, -f}};

Let us take a pseudo-polar grid with 1 ≤ r ≤ 2, \(TraditionalForm\`\(-π\) ≤ θ ≤ π\):

grid = {{1, 2}, {-Pi, Pi}};

A plot of the Poincaré map at time 4:

PoincareMapPlot[field, {t, 0, 4}, grid, GridPoints -> {7, 25}, 
        GridParameterization -> myParameterization, 
        AspectRatio -> Automatic];

If  ϕ(x, y) the Poincaré map at time 4, consider the real function d(x, 
          y) := ∥ϕ(x, y) - (x, 
              y) (∥ ^2). A very small value for d means a solution pretty much periodic with period 4. Here is a 3D plot of  d(x, y):

PoincareMapFixedPointPlot3D[field, {t, 0, 4}, grid, 
        GridPoints -> {10, 41}, PlotPoints -> {20, 41}, 
        GridParameterization -> myParameterization, 
        AxesLabel -> {"x", "y", "z"}, 
        PlotRange -> {0, 1}];

With a density plot of Min[d(x, y), ϵ] the regions where d(x, y)is smaller than ϵ will stand out as grey or black against a white background. In the following picture ϵ = 1/100. The parameters on the two axes are not x, y but r, θ:

PoincareMapFixedPointDensityPlot[field, {t, 0, 4}, grid, 1/100, 
        GridPoints -> {20, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> myParameterization, 
        FrameLabel -> {"r", "θ"}];

We make out a grey-black stripe that looks like the graph of a set  r = constant. Let us zoom in on the point r = 1.72, θ = 0:

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 4}, {{1.6, 1.8}, {- .2,  .2}}, 1/1000, 
        GridPoints -> {20, 21}, PlotPoints -> {20, 20}, 
        GridParameterization -> myParameterization, 
        FrameLabel -> {"r", "θ"}];

Let us minimize d near the point  r = 1.72, θ = 0 with the standard minimization routines:

{squareDistance, initialData} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 4}, {{1.6, 1.8}, {- .2,  .2}}, {r, 1.72, 
          1.73}, {θ, - .01,  .01}, GridPoints -> {40, 21}, 
        GridParameterization -> myParameterization]

Here is a plot of the solution corresponding the fixed point found above, seen as a function t ↦ x(t):

{{solutionX}, graphOfX} = 
        NDSolvePlot[{x''[t] == (-f /. x -> x[t]), 
            x[0] == (myParameterization[r, θ] /. 
                    initialData)[[1]], x'[
                0] == (myParameterization[r, θ] /. 
                    initialData)[[2]]}, x, {t, 0, 4}];

A check to see how much the values at t = 0 and t = 4 are close to each other:

{x[4] - x[0], x'[4] - x'[0]} /. solutionX

A graph of our 4-periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

{{solutionXY}, phaseSpacePlot} = 
        NDSolvePlot[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == (myParameterization[r, θ] /. 
                    initialData)[[1]], 
            y[0] == (myParameterization[r, θ] /. 
                    initialData)[[2]]}, {x, y}, {t, 0, 4}, 
          AspectRatio -> Automatic, 
          HowManyArrows -> 3];