Off[General::spell, General::spell1, NDSolve::ndsz, FindMinimum::fmgz]
Off[Function::fpct];  (* because of bug in 2.2 FilterOptions *)

PlotStyle::usage = "PlotStyle is an option for Plot, ParametricPlot,
ListPlot, VisualDSolve, PhasePlot, and SystemSolutionPlot.
PlotStyle -> style specifies that all lines or points are to be generated
with the specified graphics directive, or list of graphics directives.
PlotStyle -> {{style1}, {style2}, ... } specifies that successive lines
generated should use graphics directives style1, style2, ... .";

Hue::usage = "Hue[h] is a graphics directive that specifies that
graphical objects that follow are to be displayed, if possible, in a
color corresponding to hue h. Hue[h, s, b] specifies colors in terms of
hue, saturation and brightness. Hue may also be used as a setting to
FieldColor, an option to PhasePlot.";

GrayLevel::usage = "GrayLevel[level] is a graphics directive that 
specifies the gray-level intensity with which graphical objects that
follow should be displayed. GrayLevel may also be used as a setting
to FieldColor, an option to PhasePlot.";

Method::usage = "Method is an option to Solve, related functions, and various
numerical functions, which specifies what algorithm to use
in evaluating the result. When used with VisualDSolve and PhasePlot
it can force certain elementary methods to be used. Choices are:
Euler and RungeKutta4. One cannot use these specific methods if
StayInWindow is set to True.";

BeginPackage["VisualDSolve`", 
   {"Graphics`Colors`", "Utilities`FilterOptions`"}];


PhasePlot3DVersion::usage = "for internal use only";

AuxiliaryVariables::usage = "AuxiliaryVariables is an option to
ToSystem that specifies which auxiliary variables to use. It can be a
single symbol (for single equations) or a list (for systems).";


CheckOptions::usage = "CheckOptions[sym, (opts)] checks to see if the
options in opts are legitimate for sym, or the symbols in the list sym.
If they are not, and sym is a list, a message referring to the first
faulty case is returned.";

Color::usage = "Color is an option to ProjectionPlot3D that, when
False, causes a black and white image to appear.";

ColorParametricPlot::usage = "ColorParametricPlot[f, {t, a, b}, (opts)]
generates a parametric plot whose color changes continuously as a
function of t. Options for ParametricPlot may be used. The
Speed option allows the coloring to vary with the speed of
motion instead.";

ColorParametricPlot3D::usage = "ColorParametricPlot3D[f, {t, a, b},
(opts)] shows a 3-dimensional space curve using Hue to color the
curves. The PlotStyle and PlotPoints options may be used.";

ComputeWindow::usage = "ComputeWindow is an option to VisualDSolve and
PhasePlot that works when StayInWindow is True. If set to Automatic,
then computation shuts off at the viewing window, but a setting of
{x0, x1} (VisualDSolve) or {{x0, x1}, {y0, y1}} (PhasePlot) causes
that setting to be used instead.";

ContourData::usage = "ContourData is an option to Equilibria that
accepts contour plot data, thus avoiding recomputation in cases where
that data has already been computed. It is only for internal use by
NullclinePlot.";

DirectionArrow::usage = "DirectionArrow is an option to PhasePlot
that, when True, adds arrows to show the direction of increasing t.";

DirectionField::usage = "DirectionField is an option to
VisualDSolve that when True, causes a field of slope lines
to appear.";

DirectionFieldPlot::usage =
"DirectionFieldPlot[f, {t, t0, t1, dt}, {x, x0, x1, dx}] gives
an array of lines showing the slopes defined by the differential
equation dx/dt = f. t0, t1, and dt give the range
and step-size for the t-values, and similarly for the x-values.
This is used by VisualDSolve and so the output is Graphics[] and
not shown. All Graphics options may be used. The form
DirectionFieldPlot[{fx, fy}, {x, x0, x1, dx}, {y, y0, y1, dy}]
also works, and is used for autonomous systems
{x' == fx(x, y), y' == fy(x, y)}. The lines will appear to have
the same length (or, for systems, have length proportional to the
vector), regardless of the aspect ratio used.";

EquilibriaPlotPoints::usage = "EquilibriaPlotPoints is an option to
Equilibria that sets the PlotPoints for the contour plot that
is used to generate seeds. To catch pairs of nearby equilibria it
may be necessary to increase this.";

Equilibria::usage = "Equilibria[funcs, {x, xmin, xmax},
{y, ymin, ymax}, (opts)] returns the list of points where the
funcs are 0. Options appropriate to FindRoot may be used. Also
a table of eigenvalues corresponding to the equilibrium points
is printed if ShowEigenvalues is True.";

EquilibriumPointStyle::usage = "EquilibriumPointStyle is an option
to PhasePlot that sets the style of the equilibrium points (which are
shown if ShowEquilibria is set to True). It can be a list or a single
graphics primitive.";

EquilibriumFindingMethod::usage = "EquilibriumFindingMethod is an
internal option to Equilibria to deal with nondifferentiable
functions.";

EquilibriaPrecisionGoal::usage = "EquilibriaPrecisionGoal is an
option to Equilibria that eliminates duplicates with 10^-prec.";

Euler::usage = "Euler is a choice for the Method option to
VisualDSolve and PhasePlot. Number of steps is set by the EulerSteps
option. One can't use special t-values; i.e., the initial values must
all be for t0 equalling the tmin passed in {t, tmin, tmax}. Euler[] can
be called separately and works in exactly the same way that NDSolve
does, except that the starting t0 in {t, t0, t1} must appear in an
initial condition as x[t0] == ***.";

EulerSteps::usage = "EulerSteps is an option to VisualDSolve and
PhasePlot that controls the number of steps if Euler's method
is selected.";

FastPlotting::usage = "FastPlotting is an option to VisualDSolve,
SystemSolutionPlot, and PhasePlot that, when True, causes the plot
to be drawn just by using line segments through the interpolating points.
This saves a lot of time over the adaptive plotting algorithm of
Plot (since the interpolating points are generally nicely adaptive by
construction); the curve will not be satisfactory in certain low-order
cases, so we have set the default to False. Also be aware that if
MaxSteps is large then the plots generated by this method might involve
a huge number of data points. This option cannot be used with VisualDSolve
if SymbolicSolution is set to True; and it works only if the default
setting of ParametricPlotFunction is used. PlotPoints is ignored if
FastPlotting is turned on.";

FieldColor::usage = "FieldColor is an option to DirectionFieldPlot,
PhasePlot, and PhaseLine that colors the flowing fish or vectors.
Choices are: a single color or one of Hue, GrayLevel, RandomHue,
or RandomGrayLevel. Hue and GrayLevel color the fish or vectors by
relative lengths using hues or grays. FieldColor is also an option
to PhaseLine, where it must be a single color.";

FieldLength::usage = "FieldLength is an option to PhasePlot and
DirectionFieldPlot that sets the length of the field lines or of the
flowing fish. The default setting is 0.8, which uses 80% of the available
space.";

FieldLogScale::usage = "FieldLogScale is an option to PhasePlot, PhaseLine,
and DirectionFieldPlot that scales the flowing fish or vectors
proportionally to the longest fish. This is useful when a relatively few
fish are several order of magnitudes longer than the fish in the interesting
regions of the plot. FishLogScale -> 10 scales all the fish to the same
length. Default is 1, which causes no change to the relative sizes. A typical
setting in a troublesome situation might be 2, 3, or 4.";

FieldMeshSize::usage = "FieldMeshSize is an option to VisualDSolve
and related functions that controls the mesh size of the direction
field. It can be a single integer or a pair of integers.";

FieldThickness::usage = "FieldThickness is an option to PhasePlot, PhaseLine,
and DirectionFieldPlot that sets the thickness of the field lines.";

FindAllCrossings::usage = "FindAllCrossings[f, {t, a, b}, (opts)] finds all 
places where f crosses the t-axis on the interval [a, b]. The Steps
option controls the resolution of the computation. No attempt is made
to find tangential roots.";

FindMinMethod::usage = "FindMinMethod is an internal function and option
value for the case that equilibria are sought for a nondifferentiable
function.";

FlowColor::usage = "FlowColor is an option to PhasePlot and 
FlowParametricPlot that sets the color of the fish shapes.";

FlowField::usage = "FlowField is an option to DirectionFieldPlot, PhasePlot,
and PhaseLine. It replaces the traditional arrows by curvy fish shapes,
randomly placed.";

FlowParametricPlot::usage = "FlowParametricPlot[f, {t, a, b},
{x, xmin, xmax}, {y, ymin, ymax}] generates a parametric plot whose
thickness increases with time to indicate direction of flow and speed of
flow.";

FlowThickness::usage = "FlowThickness is an option to PhasePlot, 
DirectionFieldPlot and PhaseLine that, when fish are used, scales the head
size, and so the overall thickness, of the fish. Default is 1.";

FreehandAttempt::usage = "FreehandAttempt[eqn, {t, tmin, tmax},
{x, xmin, xmax}, drawPts, (opts)] compares the curve determined by
drawPts to the actual solution (in red) generated from the first
points in drawPts. It works only for VisualDSolve examples, and
passes its options to VisualDSolve.";

FreeOfSetQ::usage = "FreeOfSetQ[expr] checks to see if Set (=) appears.";

GetPts::usage = "GetPts[interFcn[t]] returns the set of
y-interpolating points in the InterpolatingFunction interFcn, when
viewed as a function y of t.";

GetPtsFull::usage = "GetPtsFull[interFcn[t]] returns the set of
interpolating points (t and y) in the InterpolatingFunction interFcn, 
when viewed as a function y of t.";

GrayShading::usage = "GrayShading is an option to ColorParametricPlot that,
when set to True, causes the colors to be replaced by shades of gray.";

Grid::usage = "Grid[n] is a possible setting of the InitialValues
option for VisualDSolve and causes an n-by-n grid to be used.";

InflectionCurve::usage = "InflectionCurve is an option to VisualDSolve
that causes the inflection curves to be drawn.";

InflectionStyle::usage = "InflectionStyle is an option to
VisualDSolve that specifies the style of the inflection lines.";

InitialPointStyle::usage = "InitialPointStyle is an option to VisualDSolve
and related functions that specifies the style of the initial-value points.
It can be a list or a single graphics primitive.";

InitialValues::usage = "InitialValues is an option to VisualDSolve
and related functions that sets the initial values. To specify
one initial value, use {t0, x0} (for systems, use {t0, x0, y0, ...}).
To specify several, use {{t0, x0}, {t1, x1}, ...}. However, one can
also attach different t-ranges for efficiency (this does not work 
with SystemSolutionPlot): for one curve use {t0, x0, {tmin, tmax}}
(or, for a PhasePlot of a system, {t0, x0, y0, (z0,...), {tmin, tmax}});
for several initial values, use several lists having this form. In the
case of systems the t-value can be omitted, and it will be set to the
lower bound in the main t-iterator. The possibility of varying the t-range
makes NDSolve more efficient, since the ranges can be tweaked until they
cause solutions that exist only in the final window; but
StayInWindow->True accomplishes this too. Also, in VisualDSolve 
a setting of Grid[n] is a possibility for InitialValues."; 

IsoclinePlotPoints::usage = "IsoclinePlotPoints is an option to
VisualDSolve that sets PlotPoints for the contour functions called
for isoclines or inflection curves.";

Isoclines::usage = "Isoclines is an option to VisualDSolve that
causes positions of constant slope to be joined.";

IsoclineShading::usage = "IsoclineShading is an option to
VisualDSolve that shows a background shaded in gray according to
the slopes. In its default form there are just two shades of gray,
since the default setting for IsoclineValues is just {0}. This will
generally be good enough, though you can get more regions by
specifying Contours to be several values, or an integer, in which
case a uniform set of levels is used.";

IsoclineStyle::usage = "IsoclineStyle is an option to VisualDSolve
that specifies the style of the isoclines that are specified in
IsoclineValues.";

IsoclineValues::usage = "IsoclineValues is an option to VisualDSolve
that specifies which slopes are to have isoclines. Default is {0}. If
it is an integer, then that many equal-spaced contours are drawn.";

LastPoint::usage = "LastPoint is a variable name that is used to store
the last point(s) of a solution.";

Logarithmic::usage = "Logarithmic is a value to the PlotType option
that requests a logarithmic (base 10) plot of the error.";

MainThickness::usage = "MainThickness is an option to ProjectionPlot3D
that sets the thickness, as a percentage of image width, of the
main spacecurve.";

MinorThickness::usage = "MinorThickness is an option to ProjectionPlot3D
that sets the thickness for the projected curves, as a percentage of image
width.";

NullclinePlot::usage =
"NullclinePlot[{f1, f2}, {x, xmin, xmax}, {y, ymin, ymax}] produces,
according to the setting of its options, the graphs of f1 = 0 and f2 = 0,
a gray background that indicates the four possibilities for
{Sign[f1], Sign[f2]}, and dots on the simultaneous roots (the equilibrium
points). It works only on orbits for autonomous systems of two equations.";

NullclinePlotPoints::usage = "NullclinePlotPoints is an option to
NullclinePlot that controls the resolution of the contour lines.";

Nullclines::usage = "Nullclines is an option to NullclinePlot and
PhasePlot that causes the nullclines to be shown.";

NullclineShading::usage = "NullclineShading is an option to NullclinePlot
and PhasePlot that causes a coded gray background to show.";

NullclineShadingMethod::usage = "NullclineShadingMethod is an
option to NullclinePlot that controls the algorithm used to
shade the regions. Possibilities are Automatic or SumOfSigns.";

NullclineStyle::usage = "NullclineStyle is an option to NullclinePlot
that sets the styles of the two nullclines. Thus it must be set to
a pair, the first entry being a list of graphics primitives
for the x-nullcline, and the second a list for the y-nullcline."; 

NumberFish::usage = "NumberFish is an option to PhaseLine and 
FlowParametricPlot that controls the number of fish.";

ParametricPlotFunction::usage = "ParametricPlotFunction is an option
to PhasePlot that, if set to ColorParametricPlot, causes
a routine called ColorParametricPlot to be used to form the orbit(s).
This colors the curve as the parameter increases. Another possible setting
is FlowParametricPlot, which causes the FlowParametricPlot function to
draw the orbit.";

PhaseLine::usage = "PhaseLine[eqn, x[t], {x, xmin, xmax}, (opts)]
shows a one-dimensional phase portrait of a first-order, autonomous
differential equation.";

PhasePlot::usage = "PhasePlot[eqn, unks, {t, tmin, tmax}, {x, xmin, xmax},
{y, ymin, ymax}, (opts)] shows the x-y orbit of the system eqn,
which can be a system of 2 or more equations (not necessarily autonomous).
Any options of ParametricPlot, Graphics, or NDSolve may be used; and there
are LOTS of special options as well. Initial values must be specified by the
InitialValues option. If there are 3 or more equations and the x-y iterators
are replaced by 3 iterators, then the orbits are given as space curves, and
options appropriate for 3D graphics may be used. In the 3D case, if the
default ParametricPlotFunction is used, then PlotPoints will be set
to 100, and can be changed. The plot is simply a uniform sequence
of 100 line segments. If only a vector field or nullcline plot
is wanted, the t-iterator, which is then irrelevant, may be omitted.";

PlaneResolution::usage = "PlaneResolution is an option to ProjectionPlot3D
that sets the grid of the three projection planes.";

PlotLabels::usage = "PlotLabels is an option to SystemSolutionPlot
that is used when more than one plot is generated and accepts a
list of plotlabels, one for each plot.";

PlotRanges::usage = "PlotRanges is an option to SystemSolutionPlot
that is used when more than one plot is generated. It must be a list
of plot ranges, one for each plot.";

PlotType::usage = "PlotType is an option to ResidualPlot that
allows a Logarithmic plot to be specified.";

PlotVariables::usage = "PlotVariables is an option to SystemSolutionPlot
that specifies which variables are to be plotted.";

PoincareSection::usage = "PoincareSection[{eq1, eq2},
{t, tmin, tmax}, {x, xmin, xmax}, {y, ymin, ymax}, {x0, y0}, (opts)]
returns a Poincare section of a system of two first-order ODEs, with the
time interval determined by the TimeInterval option. {x0, y0} are the
initial values. Equations are entered using x'[t], x[t], and so on. It
also works on a single second-order equation via: PoincareSection[
eqn, {t, tmin, tmax}, {x, xmin, xmax}, {x', ymin, ymax}, {x0, xp0}].";

PointStyle::usage = "PointStyle is an option to PoincareSection that sets
the style in which the points are shown.";

ProjectionPlot3D::usage = "ProjectionPlot3D[{f, g}, {t, tmin, tmax},
{xmin, xmax}, {ymin, ymax}] generates the space curve {f, g, t}
together with its three planar projections.";

Rainbow::usage = "Rainbow is an option to VisualDSolve, FlowParametricPlot,
PhasePlot, and
SystemSolutionPlot that causes the curves corresponding to different
initial values to appear in a variety of hues. It also works with 
PoincareSection, in which case different hues are used for the points.";

RandomGrayLevel::usage = "RandomGrayLevel is a possible value
to the FieldColor option for DirectionFieldPlot and PhasePlot
which colors the fish or vectors with random shades of gray.";

RandomHue::usage = "RandomHue is a possible value
to the FieldColor option for DirectionFieldPlot and PhasePlot
which colors the fish or vectors with random hues.";

ResidualPlot::usage = "ResidualPlot[eqn, soln, unk, {t, tmin, tmax}, (opts)]
plugs a potential solution of a single (first- or second-order) ODE
into the equation, subtracts, and plots the difference. unk is an expression
representing the dependent function of t, typically x[t].";

RKSteps::usage = "RKSteps is an option to RungeKutta4 and routines
that call it (VisualDSolve, PhasePlot, SecondOrderPlot) that controls
the number of steps if the Runge-Kutta method is selected.";

RungeKutta4::usage = "RungeKutta4 is a choice of method for VisualDSolve
and PhasePlot, and forces the standard 4th-order Runge-Kutta method. Steps
are set by the RKSteps option. RungeKutta4 can be called as a separate
function; it works in exactly the same way that NDSolve does, except that
the starting t0 in {t, t0, t1} must appear in an initial condition as
x[t0] == ***.";

SaveLastPoint::usage = "SaveLastPoint is an option to VisualDSolve,
PhasePlot, and SystemSolutionPlot that saves the final point(s) in LastPoint.
";

SaveSolution::usage = "SaveSolution is an option to VisualDSolve
and related functions that stores the solutions in a variable called
'Solution' so that it can be accessed after the use of VisualDSolve.
It does not work if symbolic solutions are sought. In the case of
PhasePlot, only the plotted solution are saved; i.e., if only two
variables out of, say, 4, are plotted, only the plotted ones are
saved. The saved object is an expression in t, not a pure function;
to evaluate it at t0, use Solution /. t -> t0.";

SecondOrderPlot::usage = "SecondOrderPlot[eqn, unks, {t, tmin, tmax}, 
(opts)] produces a plot of the solution to eqn, a second-order equation or
system of such. Plotting variables are specified by the PlotVariables
option, and x' may be used. This calls SystemSolutionPlot, and so plots
may be requested in any of the forms that that routine produces. The default
is to show graphs of all the functions and their derivatives on a single
plot. If exactly two iterators such as {x, 0, 1}, {x', 0, 1} are included
before the options, then PhasePlot is called and a phase plot is shown;
options for PhasePlot may be used. The initial values must be given
for the main variables first (unks), then the implicit auxiliary variables.";

SeedsOnly::usage = "SeedsOnly is an option to Equilibria that
returns rough values of the points. It is faster and often good
enough for visualization. The main point is that this allows
the nullclines to be split up into pieces at the equilibrium points.";

Segments::usage = "Segments is an option to FlowParametricPlot that
sets the number of segments used to draw each fish.";

ShowAuxiliaryVariables::usage = "ShowAuxiliaryVariables is an option to
ToSystem that causes the auxiliary variables to be returned as a second
argument.";

ShowEigenvalues::usage = "ShowEigenvalues is an option to Equilibria
that causes a table of eigenvalues to be displayed.";

ShowInitialValues::usage = "ShowInitialValues is an option to VisualDSolve
and related functions that when set to True, causes the initial-value
points to be superimposed on the image.";

SavePoints::usage = "SavePoints is an option to PoincareSection that
causes the points to be returned.";

ShowEquilibria::usage = "ShowEquilibria is an option to
PhasePlot and NullclinePlot that causes the equilibrium
points to be shown.";

ShowNumberOfSteps::usage = "ShowNumberOfSteps is an option to VisualDSolve,
PhasePlot, and SystemSolutionPlot that yields a printed message with the 
number of steps that the numerical solver used (same as the number of 
interpolating points in the interpolating function representing the solution).
";

SolutionName::usage = "SolutionName is an option to VisualDSolve
and related functions that accepts a string (!), which will become
a global variable used to store the solution if SaveSolution is True.";

Source::usage = "Source is an option to DirectionFieldPlot that tells if
the source was a single equation or system. For internal use only.";

Speed::usage = "Speed is an option to ColorParametricPlot that,
when True, causes the coloring to run from blue (slowest) to red
(fastest) corresponding to the speed of the particle.
If GrayShading and Speed are both True, then the shading goes from
black (slowest) to light gray (fastest).";

StayInWindow::usage = "StayInWindow is an option to
PhasePlot and VisualDSolve that, when True, causes the plots to stop
when they hit the window frame (or the max value in the t-iterator,
or the setting of ComputeWindow). It also causes the orbits to expand
in both directions from the initial point(s). When this setting is
used, Method must be set to the default of Automatic. When
StayInWindow is False, then the specified t-values are used. If any
initial values come with specific t-limits, then they override a
StayInWindow request. This allows specific t-limits to be used
for periodic solutions that will never leave the window. The exit
data, in a form suitable for later use with InitialValues, is
returned as output.";

Steps::usage = "Steps is an option to FindAllCrossings that specifies
the number of times the function is evaluated.";

SumOfSigns::usage = "SumOfSigns is a value to the
NullclineShadingMethod that forces a different, simpler, but less
highly resolved, method to be sued for shading nullclines.";

SymbolicSolution::usage = "SymbolicSolution is an option to
VisualDSolve that requests the solver to precede the numerical
approach with an attempt at finding a symbolic solution.";

SystemSolutionPlot::usage = "SystemSolutionPlot[equations, unknowns,
{t, tmin, tmax}, (opts)] generates plots of x vs. t for every x that is
a plotting variable (controlled by the PlotVariables option. Initial
values must be included, and can have the forms {t0, x0, y0, ...} or
just {x0, y0, ...} or a list of such; missing t-values are assumed to be
tmin.";

Thinning::usage = "Thinning is an option to FreehandAttempt that
thins out the user's data points to every third point, to speed up
the interpolation to get the user's function.";

TickLabelsOnly::usage = "TickLabelsOnly is an option to PhasePlot
that causes the tick labels to appear, but no tick lines. This is
useful because in the case of a shaded nullcline plot we have
arranged for the tick lines to be on the outside of the axis. If
the user does not want this, then this option allows him or her to
turn it off, but retain the tick labels.";

TimeInterval::usage = "TimeInterval is an option to PoincareSection
that sets the time interval in the section. The default setting of
Automatic causes the interval to be set to 1/100 of the t-interval.";

TimeScale::usage = "TimeScale is an option to PhaseLine that sets the
time period that is used by NDSolve to calculate the flow from each
initial condition. The initial conditions are chosen internally.
Specifying a larger or smaller TimeScale affects how long each 
arrow or fish will be. The time-length of each fish will be 
TimeScale/NumberFish.";

ToSystem::usage = "ToSystem[eqn, unk, t] turns a second-order DE to
a system of two first-order DEs. If eqn is a list of second-order DEs
and unk is a list of unknown functions, then a system of 2n equations
is returned.";

VectorField::usage = "VectorField is an option to PhasePlot
(autonomous case only) that, when True, causes a field of vectors
to appear. Fish shapes are possible for phase plane plots via the
FlowField option.";

VisualDSolve::usage =
"VisualDSolve[equation, {t, tmin, tmax}, {x, xmin, xmax}, (opts)]
creates an image that shows the solutions to a first-order ODE
superimposed on a field of slope lines. There are LOTS of options!
And all of the options for Graphics, Plot, ContourPlot, and NDSolve
may be used as well. The equation must be input in the standard
form such as x'[t] + x[t] == t + x[t]^2. Note also that the first
iterator defines the independent variable and the other iterator
defines the dependent variable.";

WindowShade::usage = "WindowShade is an option to VisualDSolve,
PhasePlot, and PoincareSection that, if set to a color or gray
shade, places a shaded rectangle in the viewing window.";

VisualDSolveVersionNumber = "2.01";

CheckOptions::badopt = "`` is not a valid option to this usage of ``.";
PhasePlot::ntatn = "Direction fields can be shown only for autonomous systems.
";
PhasePlot::noinv = "You have not included any initial values. Neither a 
direction field nor a nullcline plot can be shown because the system is not 
autonomous. Please specify initial values with the InitialValues option.";
PhasePlot::noder = "FlowField could not be calculated because the function 
could not be differentiated symbolically.";
PhasePlot::nothng = "There are no initial values and no requests via the 
Nullclines, FlowField, VectorField, NullclineShading, or ShowEquilibria 
options. You must specify something to show.";
PhasePlot::incons = "You cannot include an InitialValues setting without 
specifying a dependent variable and bounds on it.";
PhasePlot::bdordr = "The dependent variables must be in the proper form and 
same order wherever they occur (equations, list of dependent variables, 
plotting window iterators).";
PoincareSection::mxst = "Maximum number of steps was reached, so 
PoincareSection cannot proceed.";

SecondOrderPlot::badvbl = "PlotVariables must be set to All or a list of 
lists. For example, a single plot of x and x' vs. t comes from the setting {{
x', x}}; separate plots come from {{x'}, {x}}.";
SystemSolutionPlot::noinit = "Initial values are required. Use the 
InitialValues option.";

SystemSolutionPlot::nopsty = "Customized plot styles are allowed only if a 
single plot is requested.";
SystemSolutionPlot::badln = "The initial conditions do not have length `` as 
they must to match the number of variables.";
SystemSolutionPlot::badvar = "The plotting variables must be subsets of the 
unknown variable names.";
SystemSolutionPlot::norain = "The Rainbow option can be used only when a 
single plot is requested.";
VDS::insarg = "`` requires at least `` arguments.";
VDS::badarg = "The arguments are not in the correct form. Check the usage 
message. Common errors are: Omitting the list of unknown functions, having 
inconsistent variable names in different parts, omitting x'[t] or the [t] part
 in the equations, using variables that have a numerical value assigned to 
them."; 
VDS::nonopt = "Options expected beyond position ``. An option must be a rule 
or a list of rules.";
VDS::nonopt1 = "Options expected after the lists. An option must be a rule or 
a list of rules.";
VDS::nonnum = "A number is expected at position ``.";
VDS::notlist = "The first argument must be a list.";
VDS::badopt = "`` is not a valid setting for ``.";
VDS::mxstep = "NDSolve ran into a MaxSteps limitation for the initial value ``
 and only computed the orbit out to `` = ``. If a larger domain is needed, 
increase the MaxSteps option setting.";
VDS::badform = "The equations must be in terms of functions of ``.";

VisualDSolve::incons = "The StayInWindow and SymbolicSolution options cannot 
both be True, because StayInWindow works only for numerical solutions.";
VisualDSolve::incons2 = "SaveLastPoint works only when SymbolicSolution is 
False.";
VisualDSolve::nothng = "There are no initial values and no requests via the 
DirectionField, Isoclines, or InflectionCurve options. You must specify 
something to show.";
VisualDSolve::incons1 = "The SolutionName (or SaveSolution) options are meant 
for numerical solutions, and so cannot be used if SymbolicSolution is True.";
VisualDSolve::badinput = "`1` has = where it should have ==.";
VisualDSolve::bdmet = "If StayInWindow is set to True you must use the default
 setting of Method.";
VisualDSolve::badname ="SolutionName must be a string (i.e., name within 
quotes).";
VisualDSolve::failure = "NDSolve failed";
VisualDSolve::manycs = "The symbolic solution has more than one case, so 
numerical techniques are being used for the plotting.";
VisualDSolve::nosym = "The equation could not be solved symbolically.";

Begin["`Private`"];

SetAttributes[{VisualDSolve, PhasePlot, FreeOfSetQ, 
 SecondOrderPlot, PhaseLine, SystemSolutionPlot}, HoldAll];
 
ver3Q = ($VersionNumber >= 3)

(* fixes the style to a list of lists *)

fixSty[s_, n_] := Module[{solSty = s /. Automatic -> {}},
If[!ListQ[solSty], solSty = {solSty}];
If[solSty == {} || Or @@ (!ListQ[#] & /@ solSty), 
  Array[solSty &, n], solSty]]
    
(* isolates the derivatives *)
NormalForm[eqn_, unk_] :=
  Solve[eqn, D[unk, unk[[-1]]]][[1,1]] /. Rule -> Equal

NormalForm[eqns_List, unks_] :=
  First@Solve[eqns, D[unks, unks[[1,-1]]]] /. Rule -> Equal

(* this allows for ticks outside the axes; it uses a vble internal to
PhasePlot   *)
fixTicks[tl_] := (fac = If[ticLO === True, 0, 1];
Map[{#, If[# == Round[#], Round[#], #], {0, fac 0.01}}&, tl, {2}]);

g = Compile[{a, b, c, d}, If[d==0 && c==0, 10^(7),
         a b / Sqrt[(b c)^2 + (a d)^2]]];

gp = Compile[{a, b, c, d}, If[d==0 && c==0, 10^(7),
         Sqrt[a b / Sqrt[(b c)^2 + (a d)^2]]/2]];
  
h[x_] := Which[
   x >= 1   - .2,  .32, 
   x >= 2/3 - .1,  .6,
   x >= 1/2 - .1,  1,
   x >= 1/3 - .17, .8, 
   True,           .46]
     
elimDups1D[data_, tol_] := Module[{d = Sort[data]},
 Union[FoldList[If[#2 < tol + #1, #1, #2]&, First[d], Rest[d]]]]

elimDups2DOrdered[data_, tol_] := 
( FoldList[If[(#2 - #1).(#2 - #1) < tol^2, #1, #2]&,
  First[data], Rest[data]] //. 
    {aa___, {xx_?NumberQ, yy_}, {xx_, yy_}, cc___} :> {aa, {xx, yy}, cc})
    
normsq[x_, set_] := Min[(x-#).(x-#) & /@ set]

elimDups2D[data_, tol_] := Fold[If[normsq[#2, #1] > tol^2,
  Append[#1,#2], #1] &, {First[data]}, Rest[data]]

CheckOptions[sym_Symbol, opts___Rule] := CheckOptions[{sym}, opts]

CheckOptions[sym_List, opts___Rule] := 
 Check[Scan[If[!MemberQ[(First /@ (Union @@ (Options /@ sym))), #],
  Message[CheckOptions::badopt, #, First[sym]]]&, 
    First /@ {opts}], False]

(* GetPts gets range data for interpolating fcns,
   GetPtsFull gets domain and range data.
   getDomain gets the domain of an Interpolating Function,
   getDomPts gets the domain points of IF[t]   *)

Clear[GetPts, GetPtsFull];

getDomain[Ifcn_[_]] := Flatten[Ifcn[[1]]]
getDomain[Ifcn_] := Flatten[Ifcn[[1]]]

getDomPts[fcn_[_]] := First /@ fcn[[2]] /;  !ver3Q
getDomPts[fcn_[_]] := fcn[[3, 1]] /;  ver3Q

GetPts[fcn_[_]] :=  #[[3, 1]] & /@ fcn[[2]] /; !ver3Q

GetPts[fcn_[tt_], {ti_, tf_}] := 
   Last /@ GetPtsFull[fcn[tt], {ti, tf}]


GetPtsFull[fcn_[_]] :=  {#[[1]], #[[3, 1]]} & /@ fcn[[2]] /; !ver3Q

GetPtsFull[fcn_[tt_], {ti_, tf_}] := Select[GetPtsFull[
  fcn[tt]], ti <= #[[1]] <= tf &]

					
GetPts[fcn_[_]] :=  First /@ fcn[[4]]/; ver3Q
GetPtsFull[fcn_[_]] :=   Transpose[{fcn[[3, 1]], First /@ fcn[[4]]}] /; ver3Q


(* ****************** FreeOfSetQ  *************  *)

FreeOfSetQ[equation_] := 
   If[FreeQ[Hold[equation], Set], True,
    Message[VisualDSolve::badinput, HoldForm[equation]]; False];


  (*   ***********  ColorParametricPlot  ***********  *)

Options[ColorParametricPlot] = Join[
 {GrayShading->False, Speed->False, PlotStyle->{AbsoluteThickness[1]}},
  Options[ParametricPlot]];

ColorParametricPlot[f_, {t_, a_, b_},  opts___] := Module[{i = 0, g},
  {gr, styles, sp} = {GrayShading, PlotStyle, Speed} /. {opts} /. 
    Options[ColorParametricPlot];

 data = Table[Null, {500}];
 options = Join[{FilterOptions[ParametricPlot, opts]}, 
  {Compiled->False, DisplayFunction->Identity}];
 options1 = Join[{FilterOptions[Graphics, opts]},
     {PlotRange->All, AspectRatio->Automatic}];
 norm[v_] := N[Sqrt[v.v]];
 
 If[sp, speed[tt_] := Evaluate[norm[D[f, t] /. t->tt]];
    {speedMin, speedMax} = (ll = Last[Transpose[Plot[speed[tt], {tt, a, b}, 
               PlotRange->All, DisplayFunction->Identity][[1,1,1,1]]]]; 
   {Min[ll], Max[ll]})];

g[s_] := (fval = N[f /. t->s];
           If[i == Length[data], data = Join[data, Array[Null &, 500]]];
           data[[++i]] = {s, fval}; fval);  

ParametricPlot[g[t], {t,a,b}, DisplayFunction->Identity, Evaluate[options]];
data = Partition[Sort[Take[data, i]], 2, 1];  
 
{c1, c2, func1, func2, aa, bb, sFunc, d} = Which[
    gr && !sp, 
    {0.8, N[Pi]/2, GrayLevel[Min[Max[#, 0], 1]]&, Sin, a, b, Identity, 0}, 
   !gr && !sp, 
   	{0.6,   1,  Hue,  	(1-#)&, a,   b,  Identity, 0},
   !gr &&  sp,
   	{0.5, 1,     Hue,  	Identity, speedMin, speedMax, speed, 0.5},
    gr && sp,
    {0.8, N[Pi/2], GrayLevel[Min[Max[#,0], 1]]&, Sin, speedMin, 
        speedMax, speed, 0}];

 Show[Graphics[{Sequence@@Flatten[{styles}], data /. {{s_, pt_List}, next_} :>
 
       {func1[d + c1 * func2[N[c2 * (sFunc[s] - aa)/(bb - aa)]]], 
            Line[{pt, next[[2]]}]}}],  Evaluate[options1], 
            Axes->Automatic, PlotRange->All  ]]


  
     (*  ***********  ColorParametricPlot3D  ********  *)

Options[ColorParametricPlot3D] = Append[Options[ParametricPlot3D], 
  PlotStyle->{}];

ColorParametricPlot3D[f_List, {t_, tmin_, tmax_}, opts___] :=
Module[{ffx, ffy, ffz},
  ffx = Compile[{t}, Evaluate[f[[1]]]];
  ffy = Compile[{t}, Evaluate[f[[2]]]]; 
  ffz = Compile[{t}, Evaluate[f[[3]]]]; 
  
{pSty, tSteps} = {PlotStyle, PlotPoints} /. {opts} /. 
    Options[ColorParametricPlot3D];
tSteps = tSteps /. Automatic->75;
If[Head[pSty] =!= List, pSty = {pSty}];
pSty = pSty /. {{seq___}} :> {seq};
(* preceding because PhaePlot might use too many {} *)

Show[Graphics3D[Append[pSty, MapIndexed[ 
 {Hue[0.6 (1 - #2[[1]]/tSteps)],  Line[#1]} &, 
  Partition[ {ffx[#], ffy[#], ffz[#]} & /@
      N[Range[tmin, tmax, (tmax - tmin)/(tSteps - 1)]], 2, 1]]],
Evaluate[FilterOptions[Graphics3D, opts]],
DisplayFunction->Identity]]]

  

Off[CompiledFunction::cfn, CompiledFunction::cfr,
   Power::infy, Infinity::indet];


                (* ProjectionPlot3D *)

Options[ProjectionPlot3D] = {MainThickness->0.008, 
MinorThickness->0.001,
PlotPoints->40, PlaneResolution->7, Color->True};

ProjectionPlot3D[{f_, g_}, {t_, tmin_, tmax_}, 
  {xmin_, xmax_}, {ymin_, ymax_}, opts___] := (

{minorth, th, pp, pres, colQ} = {MinorThickness,
  MainThickness, PlotPoints, PlaneResolution, Color} /. 
    {opts} /. Options[ProjectionPlot3D];


colors = If[colQ, 
 {RGBColor[0,1,1], RGBColor[.5,1,0.83], RGBColor[.2, .36, 1], 
  RGBColor[1,0,0], RGBColor[1,0,1], RGBColor[1, 1, 0]},
    {GrayLevel[1], GrayLevel[1], GrayLevel[1], GrayLevel[.5],
      GrayLevel[.5],  GrayLevel[.5]}];

Show[
 ParametricPlot3D[{x, y, tmin, colors[[1]]},
   {x, xmin, xmax}, {y, ymin, ymax}, PlotPoints->pres, 
   Lighting->False, DisplayFunction->Identity], 

 ParametricPlot3D[{xmin, y, t, colors[[2]]},
   {t, tmin, tmax}, {y, ymin, ymax}, PlotPoints->pres, 
			Lighting->False, DisplayFunction->Identity], 
 
 ParametricPlot3D[{x, ymax, t, colors[[3]]},
   {t, tmin, tmax}, {x, xmin, xmax}, PlotPoints->pres, 
			Lighting->False, DisplayFunction->Identity], 

 Graphics3D[{
 {Thickness[th],
   Line[Append[
     Table[{f, g, t}, 
       {t, N@tmin, N[tmax - (tmax - tmin)/pp], N[(tmax-tmin)/pp]}],
       {f, g, t} /. t->tmax]]},
  
 Thickness[minorth], {colors[[4]], Line[{#[[1]], #[[2]], tmin} & /@  
   ParametricPlot[{f, g}, {t, tmin, tmax}, 
     DisplayFunction->Identity][[1,1,1,1]]] },
     
 {colors[[6]], Line[{#[[2]], ymax, #[[1]]} & /@  
   Plot[f, {t, tmin, tmax}, 
     DisplayFunction->Identity][[1,1,1,1]]] },
     
  {colors[[5]], Line[{xmin, #[[2]], #[[1]]} & /@  
   Plot[g, {t, tmin, tmax}, 
     DisplayFunction->Identity][[1,1,1,1]]] }  }],
  FilterOptions[Graphics3D, opts],
  Boxed->False,
  AxesLabel->{None, None, t}])


(* ****************** DirectionFieldPlot  *************  *)

Options[DirectionFieldPlot] = {
  Compiled -> False,
  FieldColor->GrayLevel[0],
  FieldLength->.8,
  FieldThickness->Thickness[0.004],
  FlowField->False,
  AspectRatio->1/GoldenRatio,
  PlotRange->Automatic,
  FlowThickness->1.,
  FieldLogScale->1.,
  Source->PhasePlot};

DirectionFieldPlot[f_ /; Head[f] =!= List, a__] := 
            DirectionFieldPlot[{1, f}, a, Source -> FirstOrder]

DirectionFieldPlot[f_List, {x_, tx0_, tx1_, tdx_}, {y_, ty0_, ty1_, tdy_}, 
opts___] :=
 Module[{fx, fy},
  {col, length, thick, asp, plotR, useFish,
    fishl, fhead, sou, comQ} = 
  {FieldColor, FieldLength, FieldThickness, AspectRatio, 
  PlotRange, FlowField, FieldLogScale, 
  FlowThickness, Source, Compiled} /. {opts} /. 
     Options[DirectionFieldPlot];

{x0,x1,dx,y0,y1,dy} = N[{tx0,tx1,tdx,ty0,ty1,tdy}];           
{fx, fy} = If[comQ, Compile, Function][{x, y}, #] & /@ f;
fnew = f /. {Abs[zz_] :> If[zz<0,-1,1] zz,
            Sign[zz_] :> If[zz<0,-1,1]};
{fxp, fyp} = If[comQ, Compile, Function][{x, y}, #] & /@ (Outer[D, fnew, {x,
y}].fnew);

If[!FreeQ[{fxp, fyp}, Derivative], Message[PhasePlot::noder]];

useFish = useFish && sou === PhasePlot && FreeQ[{fxp,fyp},Derivative];
showarrows = (sou === PhasePlot);

{{xx0, xx1}, {yy0, yy1}} = 
      N[plotR /. Automatic->{{x0, x1}, {y0, y1}}] +
      N[{{dx/4,-dx/4},{dy/4,-dy/4}}];
asp = N[asp /. Automatic->(y1-y0)/(x1-x0)];
 
{bx, by} =  N[{xx1 - xx0, (yy1-yy0)} * If[asp<1, {asp, 1}, 
   {1, 1/asp} ]/(mm=Max[(xx1-xx0)/dx, (yy1-yy0)/dy])];

unitSlope[xx_Real, yy_] := ( 
 {tx, ty} =  {fx[xx, yy], fy[xx, yy] /. 
     {Indeterminate->10^7 dy, ComplexInfinity->10^7 dy}};
    Line[{{xx, yy} - (u = (length*{tx, ty} * g[bx, by, tx, ty])/2),
           {xx, yy} + u}]);
                   
Which[col === GrayLevel,     fcol = (GrayLevel[0.8*(1-#)]&),
      col === Hue,    fcol = (Hue[(#+1.)/2.]&),
      col === RandomGrayLevel, fcol = (GrayLevel[Random[]]&),
      col === RandomHue, fcol = (Hue[Random[]]&),
      Head[col] === Function, fcol = col,
      True,        fcol = col &];
      
numfish=Round[N[(x1-x0)(y1-y0)/dx/dy]];     
fishnum=0;

If[useFish, 

perpUnitSlope[xx_Real, yy_, tx_, ty_, gg_] := 
(u = fhead*{tx, ty} * gg/8.;
pu = fhead*{-ty bx/by, tx by/bx}*gg/8.; 
txp = fxp[xx,yy]/2;
typ = fyp[xx,yy]/2;
ggp=Min[gg, gp[bx,by,txp,typ]];
flog = fishs*Exp[Log[ggp]*fishl]*length; 
{xx2,yy2} = {xx,yy} - {tx,ty}*flog + flog^2*{txp,typ};
{xx3,yy3} = {xx,yy} - {tx,ty}*2flog + 4flog^2*{txp,typ};
{xx4,yy4} = {xx,yy} - {tx,ty}*3flog + 9flog^2*{txp,typ};
If[xx4<=xx1 && xx4>=xx0 && yy4<=yy1 && yy4>=yy0,
   {fcol[N[(fishnum++)/(numfish+1)]], 
   Polygon[{ {xx,yy} + pu, {xx2,yy2} + 2pu/3 -.745u, 
     {xx3,yy3} + pu/3 -.943u, {xx4,yy4} - u,{xx3,yy3} - pu/3 -.943u, 
     {xx2,yy2} - 2pu/3 -.745u, {xx,yy} - pu}],
   Disk[{xx,yy},{bx,by}*fhead/8, {0, 2 Pi}]}] ), 
                
perpUnitSlope[xx_Real, yy_, tx_, ty_, gg_] := (
    u = length*{tx,ty}*fishs*Exp[Log[gg]*fishl];
    pu = length*{-ty bx/by, tx by/bx}*fishs*Exp[Log[gg]*fishl]/4;
    If[xx+u[[1]]<=xx1 && xx+u[[1]]>=xx0 &&  
                 yy+u[[2]]<=yy1 && yy+u[[2]]>=yy0,
      {fcol[N[(fishnum++)/(numfish+1)]],
        thick,Line[({xx, yy} + u + #) & /@
        (fhead*{pu  - .25 u , {0,0}, -pu - .25 u})],
        Line[{{xx, yy}, {xx, yy} + u}]}])];

If[useFish, 

(points = Sort[Flatten[Table[
       {aa = N[i] + Random[]*dx, bb = N[j] + Random[]*dy,
        fxaa=fx[aa,bb],fyaa=fy[aa,bb],
        g[bx,by,fxaa,fyaa]},
        {i,x0,x1-dx,dx},{j,y0,y1-dy,dy}],1],
                     (#2[[5]] < #1[[5]])&]; 
minfactor = Last[points][[5]];
midfactor = points[[Round[numfish/4]]][[5]];
fishl = fishl /. Automatic->N[Exp[Log[midfactor/minfactor/4]/
              Log[10,midfactor/minfactor]]];
fishl = N[Log[10,Max[Min[fishl,10.],1.]]];
fishs = minfactor^(1-fishl)),

(points = Join[
  Distribute[ List [
    Range[N[(x0+dx/2) /. 0->0.], (x1-dx/2), N[dx]], 
    Range[N[(y0+dy/2) /. 0->0.], (y1-dy/2), N[2dy]] ], List],
  Distribute[  List[
    Range[N[(x0+dx) /. 0->0.], (x1-dx/2), N[dx]],
    Range[N[(y0+3dy/2) /. 0->0.], (y1-dy/2), N[2 dy]]], List] ]; 
If[showarrows,
  (points = Map[ {aa=#[[1]],bb=#[[2]],
               fxaa=fx[aa,bb],fyaa=fy[aa,bb],
               g[bx,by,fxaa,fyaa]}&,points]; 
  points = Sort[points, (#2[[5]] < #1[[5]])&];     
  minfactor = Last[points][[5]];
  midfactor = points[[Round[Length[points]/4]]][[5]];
  fishl = fishl /. Automatic->N[Exp[Log[midfactor/minfactor/4]/
              Log[10,midfactor/minfactor]]];
  fishl = N[Log[10, Max[Min[fishl,10.], 1.]]];
  fishs = minfactor^(1-fishl))])];

Graphics[Join[If[!showarrows, {col, thick},{}], 
  DeleteCases[Apply[If[showarrows, perpUnitSlope, unitSlope], 
      points, 2  (* try:  {1}  *) ], Null]], 
 Evaluate[FilterOptions[Graphics, opts]]]]


Options[FindAllCrossings] = {Steps -> 113, WorkingPrecision ->
 $MachinePrecision};

FindAllCrossings[x___] :=  (Message[VDS::insarg, FindAllCrossings, 2]; 
    HoldForm[FindAllCrossings[x]]) /; Length[{x}] < 2
    
FindAllCrossings[x__] :=  (Message[VDS::nonopt, 2]; 
     HoldForm[FindAllCrossings[x]]) /; 
     Length[{x}] != 2 && Union[Head /@ Drop[{x}, 2]] =!= {Rule}

FindAllCrossings[x__] :=  (Message[VDS::badarg]; 
            HoldForm[FindAllCrossings[x]]) /; Head[{x}[[2]]] =!= List
            
FindAllCrossings[x__] :=  (Message[VDS::nonnum, {2,2}]; 
           HoldForm[FindAllCrossings[x]]) /; !NumberQ[N[{x}[[2,2]]]]
           
FindAllCrossings[x__] :=  (Message[VDS::nonnum, {2,3}]; 
           HoldForm[FindAllCrossings[x]]) /; !NumberQ[N[{x}[[2,3]]]]

FindAllCrossings[f_, {t_, a_, b_}, opts___] := Module[{},
 {wp, res} = {WorkingPrecision, Steps} /. {opts} /. Options[FindAllCrossings];


If[!IntegerQ[res], Message[VDS::badopt, res, FindAllCrossings];
      Return[HoldForm[FindAllCrossings[f, {t, a, b}, opts]]]];
fC = Compile[t, f];
s1 = Sign[fC[#] & /@  (r = Range[N@a, b, N[(b-a)/res]])];
If[!MemberQ[Abs[s1], 1], Return[{}]];
s = Drop[s1 * RotateLeft[s1], -1];
If[MemberQ[s, -1] || MemberQ[Drop[Rest[s], -1], 0], 
  N[Union[Round[1000 * 
  Join[r[[#]] & /@ First /@ Position[s1, 0], 
  Select[t /. (FindRoot[f, {t, {r[[#]], r[[#+1]]}}, 
                WorkingPrecision->wp] & /@
 First /@ Position[s, -1]), a <= # <= b &]]] / 1000]], {}]]

                     (* ***** Equilibria ******** *)

Options[Equilibria] = {ContourData -> Automatic,
  EquilibriaPlotPoints->49,
  EquilibriaPrecisionGoal->5, 
  SeedsOnly->False,
  ShowEigenvalues->False,
  EquilibriumFindingMethod->FindRoot};


Equilibria[x___] :=  (Message[VDS::insarg, Equilibria, 3]; 
    HoldForm[Equilibria[x]]) /; Length[{x}] < 3

Equilibria[x__] :=  (Message[VDS::nonopt, 3]; HoldForm[Equilibria[x]]) /; 
     Length[{x}] != 3 && Union[Head /@ Drop[{x}, 3]] =!= {Rule}

Equilibria[x__] :=  (Message[VDS::badarg]; HoldForm[Equilibria[x]]) /; 
   Head[{x}[[1]]] =!= List || Head[{x}[[2]]] =!= List || Head[{x}[[3]]] =!= 
List

Equilibria[x__] :=  (Message[VDS::nonnum, {2,2}]; 
           HoldForm[Equilibria[x]]) /; !NumberQ[N[{x}[[2,2]]]]

Equilibria[x__] :=  (Message[VDS::nonnum, {2,3}]; 
           HoldForm[Equilibria[x]]) /; !NumberQ[N[{x}[[2,3]]]]

Equilibria[x__] :=  (Message[VDS::nonnum, {3,2}]; 
           HoldForm[Equilibria[x]]) /; !NumberQ[N[{x}[[3,2]]]]

Equilibria[x__] :=  (Message[VDS::nonnum, {3,3}]; 
           HoldForm[Equilibria[x]]) /; !NumberQ[N[{x}[[3,3]]]]

(* Equilibria solves f = g = 0 as follows. Given a seed value,
it uses Newton's method via FindRoot. To generate the seeds it
uses the data from a plot of f == 0 (obtained via ContourPlot!).
It then looks at the values of g on this data set, finds the
sign changes, and uses the resulting info to generate the seeds. *)

(* Now we generate seeds, and use the preceding routine to
speed into the actual solutions. We allow contour data to be
fed in, so that when NullclinePlot calls this it does not
recompute contour data. *)

FindMinMethod[{eqn1_, eqn2_}, {x_, x0_}, {y_, y0_}, {xR_, yR_}] :=
  FindMinimum[Evaluate[eqn1^2 + eqn2^2], {x, {x0, x0 + xR/100}}, 
 {y, {y0, y0 + yR/100}}][[2]]

EigenData[funcs_List, vbles_List, ePt_] :=
Eigenvalues[(*N@*)Outer[D, funcs, vbles] /.
   Thread[Rule[vbles, ePt]]]

Equilibria[{0, 0} | {0., 0.}, __] := {} 

Equilibria[funcs_List, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, 
  opts___] := Module[
 {done = False, method, xR = xmax - xmin, yR = ymax - ymin, 
   fy = Compile[{x, y}, Evaluate@funcs[[2]]]},

 {prec, sloppy, saveEq, eigQ, method, conData} = 
    {EquilibriaPrecisionGoal, SeedsOnly, SaveEquilibria,
      ShowEigenvalues, EquilibriumFindingMethod, ContourData} /. 
        {opts} /. Options[Equilibria];

(* try Solve if funcs are polys *)
If[And @@ PolynomialQ[funcs, {x, y}] && FreeQ[funcs, _Real],
  done = (
  fail =!= (s=Check[Solve[Thread[Equal[funcs, {0,0}]],{x,y}], fail]));
  If[done = done && FreeQ[s, ToRules], eqPts = Select[
    {x, y} /. s,  NumberQ[#[[1]]]&& NumberQ[#[[2]]] &&
     xmin < N@#[[1]] < xmax &&  ymin < N@#[[2]] < ymax &];
      If[eqPts == {}, done = False]]];

If[!done && conData === Automatic, conData = First /@ Cases[Graphics[
      ContourPlot[funcs[[1]], {x, xmin, xmax}, {y, ymin, ymax},
      Contours->{0}, ContourShading->False,
      DisplayFunction->Identity,
      PlotPoints->
       (EquilibriaPlotPoints /. {opts} /. Options[Equilibria])]],
        _Line, Infinity] ];
If[!done, seeds = Flatten[ Map[#[[(1+Flatten[
       Position[
         Rest[tt = Sign[Apply[fy, #, {1}]]] * Rest[RotateRight[tt]],
          -1]])]]&, conData], 1]];         
If[!done && (sloppy || seeds == {}), eqPts = seeds; Return[eqPts]]; 
If[!done, eqPts = Chop@Select[elimDups2D[Map[
    {x, y} /. If[method === FindMinMethod,
      FindMinMethod[{funcs[[1]], funcs[[2]]},
       {x, #[[1]]}, {y, #[[2]]}, {xmax - xmin, ymax - ymin},
       Evaluate@FilterOptions[FindRoot, opts]],
     FindRoot[{funcs[[1]] == 0, funcs[[2]] == 0},
        {x, #[[1]]}, {y, #[[2]]},Evaluate@FilterOptions[FindRoot, opts]]]&,
    seeds], 10^(3-prec)],
    N@xmin < #[[1]] < N@xmax  && N@ymin < #[[2]] < N@ymax &]];
    
If[eigQ && eqPts != {}, Print[TableForm[ Chop@Transpose[
         {eqPts, EigenData[funcs, {x, y}, #]& /@ eqPts}] /.
             xx_Real :> N[xx, 2], 
         TableDepth->2, 
         TableSpacing->{1, 5},
         TableHeadings ->{None,  {"Equilibria", "Eigenvalues\n"}}]]];
   eqPts]


(*  ********   Euler  ******  *)

Options[Euler] = {WorkingPrecision->$MachinePrecision, EulerSteps->20};

Euler[eqn_, unknowns_, {t_, t0_, t1_}, opts___] := 
  Module[{g, h},
  {wp, steps} = {WorkingPrecision, EulerSteps} /. {opts} /. Options[Euler];
   fixZero = 0 -> SetAccuracy[0, wp];
   exprs = {unknowns /. y_[t] :> D[y[t], t], unknowns /. t -> (t0 /. fixZero)
 } /. 
       Flatten[ ToRules /@ N[eqn, wp] /. fixZero];
   FF = First[exprs];
  h = (t1-t0)/steps;
  g[{u_, v_}] := N[{u, v} + 
        h {1, (FF /. Flatten[Append[{Thread[unknowns->v]}, t->u]])}, wp];
  data = N[NestList[g, {t0, Last[exprs]}, steps], wp];
  data[[-1, 1]] = t1;
  If[Length[unknowns] == 1,  {{unknowns->Interpolation[data][t]}},
    {tdata, vdata} = Thread[data];
    {Thread[Rule[unknowns, 
       Table[Interpolation[Thread[{tdata, Thread[vdata][[i]]}]][t], 
     {i, Length[unknowns]}]]]}  ]]
  


(*  ********   RungeKutta4  ******  *)

Options[RungeKutta4] = {WorkingPrecision->$MachinePrecision, RKSteps->20};

RungeKutta4[eqn_, unknowns_, {t_, t0_, t1_}, opts___] := 
  Module[{g, h, FF, f1, f2, f3},

  {wp, steps} = {WorkingPrecision, RKSteps} /. {opts} /. Options[RungeKutta4];

  fixZero = 0 -> SetAccuracy[0, wp];
  exprs ={unknowns /. y_[t] :> D[y[t], t], unknowns /. 
      t->(t0 /. fixZero) }   /. 
       Flatten[ ToRules /@ N[eqn, wp] /. fixZero];
  
  h = (t1-t0)/steps;
 
  FF[u_, v_] := First[exprs] /. Flatten[Append[{Thread[unknowns->v]}, t->u]];
 
  g[{u_, v_}] := (tmid = u + h/2;
   f1 = N[h FF[u, v], wp];
   f2 = h FF[tmid, v + f1/2];
   f3 = h FF[tmid, v + f2/2];
   N[{u, v} + {h, (f1 + 2 (f2 + f3) + h FF[u + h, v + f3]) / 6}, wp]);
   
  data = N[NestList[g, {t0, Last[exprs]}, steps], wp];
  data[[-1, 1]] = t1; 
  If[Length[unknowns] == 1,  {{unknowns->Interpolation[data][t]}},
    {tdata, vdata} = Thread[data];
    {Thread[Rule[unknowns, 
       Table[Interpolation[Thread[{tdata, Thread[vdata][[i]]}]][t] , 
          {i, Length[unknowns]}]]]} ]]


Options[FlowParametricPlot] = {
  NumberFish->5,
  Segments->10,
  FlowThickness->1,
  AspectRatio->Automatic,
  FlowColor->GrayLevel[0]};

FlowParametricPlot[f_, {t_, tti_, ttf_},
        {x_, x0_, x1_}, {y_, y0_, y1_}, opts___] :=
Module[{ff},  
{col, fhead, asp, numfish, segs, rainb} = {FlowColor, FlowThickness,
   AspectRatio, NumberFish, Segments, Rainbow} /. {opts} /. Options[
FlowParametricPlot];

{ti, tf} = N[{tti, ttf}];

Which[col === GrayLevel, fcol = (GrayLevel[0.8*(1-#)]&),
      col === Hue,       fcol = (Hue[(#+1)/2]&),
      True,              fcol = col];
         
asp = N[asp /. Automatic -> (y1 - y0)/(x1 - x0)];
{bx, by} =  N[{x1 - x0, y1 - y0}*
    If[asp < 1, {asp, 1}, {1, 1/asp}]*fhead/120];

Computepu[{tx_, ty_}] := 
  {-ty bx/by, tx by/bx}*g[bx, by, tx, ty];   
  
If[Head[f[[1,0]]] === InterpolatingFunction,
   {t0, t1} = Flatten[f[[1, 0, 1]]]; 
   ff = Evaluate[{f[[1,0]][#], f[[2,0]][#]}]&,
   {t0, t1} = N[{ti, tf}];
   ff = Evaluate[{First[f] /. t->#, Last[f] /. t->#}]&];

zz = 0; zzz = 0;
points = Map[If[# < t0, zz++; t0, If[# > t1, zzz++; t1, #]]&, 
             Range[ti, tf, N[(tf - ti)/segs/numfish]]];
points = Map[ff, Drop[Drop[points, Quotient[zz, segs]*segs],
             -Quotient[zzz, segs]*segs]];
zz = Mod[zz, segs]; zzz = Mod[zzz, segs];          

pupoints =  Map[Computepu, points - RotateRight[points]];

pupoints = MapIndexed[(#1 * Mod[(#2[[1]]-1),segs]/(segs-1))&,
           Join[Table[pupoints[[zz+1]], {i, 1, zz}],
             Drop[pupoints, zz]]];
      
points = Partition[points, segs];
pupoints = Partition[pupoints, segs];

fishheads = Map[Disk[Last[#], {bx,by}, {0, 2 Pi}]&, points];
If[zzz > 1, fishheads = Drop[fishheads, -1]];

points = Apply[Plus, {Map[Join[#,Reverse[#]]&, points],
          Map[Join[#,Reverse[-1*#]]&, pupoints]}];
          
If[Head[fcol] === Function,
  points = {MapIndexed[
    {fcol[(Mod[#2[[1]]-1,segs-1]+1)/segs],Polygon[#]}&,
    Flatten[Map[MapThread[Join,{Take[#,segs - 1],
      Reverse[Take[#, 1 - segs]]}]&,
      Map[Partition[#,2,1]&, points]],1]],
      Map[{fcol[1],#}&, fishheads]}, 
  points = {fcol, Map[Polygon, points], fishheads}];
 
Graphics[points, Evaluate[FilterOptions[Graphics, opts]]]]




(* ****************** FreehandAttempt  *************  *)

Options[FreehandAttempt] = {Thinning->True};

FreehandAttempt[x___] :=  (Message[VDS::insarg, FreehandAttempt, 4]; 
    HoldForm[FreehandAttempt[x]]) /; Length[{x}] < 4
                
FreehandAttempt[x__] :=  (Message[VDS::nonnum, {r[[1,1]]+1, 2}]; 
            HoldForm[FreehandAttempt[x]]) /;
        (r=Position[(List[x][[{2,3}, 2]]), 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}

FreehandAttempt[x__] :=  (Message[VDS::nonnum, {r[[1,1]]+1, 3}]; 
            HoldForm[FreehandAttempt[x]]) /;
        (r=Position[(List[x][[{2,3}, 3]]), 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}

FreehandAttempt[eqn_Equal, {t_, tmin_, tmax_}, {x_, xmin_, xmax_}, 
   drawPts_List, opts___] := Module[{},

If[FreeQ[eqn, x'[t]], Message[VDS::badform, t]; Return[HoldForm[
FreehandAttempt[eqn, {t, tmin, tmax}, {x,xmin,xmax}, drawPts, opts]]]];

th = Thinning /. {opts} /. Options[FreehandAttempt];
d = drawPts[[Range[1, Length[drawPts], If[th, 3, 1]]]];

{fx, fy} = Interpolation[Thread[ {Range[0, 1, 1/(Length[d] - 1)], #}]] & /@
               Thread[d];
Show[VisualDSolve[eqn, {t, tmin, tmax}, {x, xmin, xmax},
      InitialValues->drawPts[[1]], opts, 
      PlotStyle->{Thickness[0.012]},
      DisplayFunction->Identity],
  ParametricPlot[{fx[tt], fy[tt]}, {tt, 0, 1},
       PlotStyle->{Dashing[{.015, .02}]},
       DisplayFunction->Identity],   
 DisplayFunction->$DisplayFunction]]


(* *******  NullclinePlot ***** *)

Options[NullclinePlot] = {
 NullclinePlotPoints->50,
 NullclineShadingMethod -> Automatic,
 ContourStyle->Automatic,
 NullclineShading->False,
 ShowEquilibria->False,
 EquilibriumPointStyle->{PointSize[0.021]},
 Nullclines->False,
 NullclineStyle->{{},{}}};

TurnPts[set_] :=  First /@ 
set[[(1+(Flatten[Position[Apply[Times,
  Partition[DeleteCases[Sign[Apply[#1 - #2 &, Partition[First /@ set, 2, 1], 
     2]], 0], 2, 1], 2  (* try {1}  *) ], -1]]))]] /; Depth[set] == 3
 
TurnPts[sets_] := Apply[Join, TurnPts /@ sets] /; Depth[sets] == 4
     
ToFunctions[set_] :=  If[#[[1,1]] > #[[-1,1]], 
   Reverse[#], #]&  /@ (set[[#]] & /@ Apply[Range,
  Partition[Prepend[Append[1+(Flatten[Position[Apply[Times,Partition[
  DeleteCases[ Sign[Apply[#1 - #2 &, 
    Partition[First /@ set, 2, 1], 2]], 0], 2,1], 2  (* try {1}  *) ], -1]]),
       Length[set]], 1], 2, 1], 2  (* try {1}  *) ]);
       
BreakAtXVal[set_, x_?NumberQ] := Module[{n = Length[set],
 i = home[First[Transpose[set]], x]},
 If[1 <= i < n,
  linInt = {{x, set[[i, 2]] + (x - set[[i,1]])*
       (set[[i+1, 2]]-set[[i,2]])/(set[[i+1, 1]] - set[[i,1]])}};
  {Join[Take[set, i], If[i<=n-1, linInt, {}]], 
   Join[If[i>0, linInt, {}], Take[set, i-n]]},  {set}] 
     ] /; Depth[set] < 4

BreakAtXVal[sets_, x_?NumberQ] := (Apply[Join, 
    BreakAtXVal[#, x] & /@ sets]) /; Depth[sets] == 4

BreakAtXVal[set_, xVals_List] := Fold[
DeleteCases[BreakAtXVal[#1, #2], {}]&, set, xVals] /; Depth[set] == 3

BreakAtXVal[sets_, xVals_List] := (Apply[Join, 
    BreakAtXVal[#, xVals] & /@  sets]) /; Depth[sets] == 4


home[xvals_, l_] :=  Infinity   /; l >= Last[xvals];
home[xvals_, l_] := -Infinity   /; l <= First[xvals];
home[{r_, s_}, _] :=    1       /; r < s;
home[xvals_, x_] := Module[{n = Length[xvals], h},
  h = Quotient[n, 2];
  Which[x < xvals[[h]],   home[Take[xvals, h], x],
        x > xvals[[h+1]], h + home[Take[xvals, h - n], x],
        True, h] ]

home1[{r_, s_}, _] := 1   /; r < s;

home1[xvals_, x_] := Module[{n = Length[xvals], h},
  h = Quotient[n, 2];
  Which[x < xvals[[h]],   home1[Take[xvals, h], x],
        x > xvals[[h+1]], h + home1[Take[xvals, h - n], x],
        True, h] ]

NullclinePlot[
  funcs_List, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, opts___] :=
 
 Module[{}, 
  {cpp, shadedQ, stabQ, cs, stabCol, showLines, shadmethod, nullcSty} =
   {NullclinePlotPoints, NullclineShading, ShowEquilibria, 
     ContourStyle, EquilibriumPointStyle, Nullclines,
       NullclineShadingMethod, NullclineStyle} /. {opts} /. 
        Options[NullclinePlot];

If[!ListQ[cs], cs = {cs}];
If[!ListQ[stabCol], stabCol = {stabCol}];
contourData =  (First /@ Cases[Graphics[
   ContourPlot[#, {x, xmin, xmax}, {y, ymin, ymax},
    Contours->{0},
    DisplayFunction->Identity, 
    Evaluate[FilterOptions[ContourPlot, opts]],
    PlotPoints->cpp,
    ContourShading->False, Frame->False]],
    _Line, Infinity]) & /@ funcs;
  
If[stabQ || (shadedQ && shadmethod === Automatic),
   stabPts = Check[Equilibria[funcs, 
     {x,xmin, xmax}, {y, ymin, ymax}, ContourData -> contourData[[1]],
         FilterOptions[Equilibria, opts],
         EquilibriumFindingMethod -> If[And @@ (FreeQ[funcs, #] & /@ {Abs, 
Sign, If, Which}),
             FindRoot, FindMinMethod]], eqFailed]];
Show[
If[shadedQ, 

If[stabPts === eqFailed || shadmethod === SumOfSigns,
DensityPlot[{2,1}.Sign[funcs], {x, xmin, xmax}, {y, ymin, ymax},
  DisplayFunction->Identity, Mesh->False,
   PlotPoints-> cpp, PlotRange->{-3, 3},
   ColorFunction->(GrayLevel[h[N[#]]]&)],

Graphics[{data3 = N[{{{xmin, ymin}, {xmax, ymin}}}];
    criticalXVals = Join[
    If[stabPts == {}, {}, First[Thread[stabPts]]],
    Union[TurnPts[jj = Join[contourData[[1]], contourData[[2]], data3]]],
    Flatten[#[[{1,-1}, 1]] & /@ jj],
    N[{xmin, xmax}]];   

criticalXVals = N[elimDups1D[criticalXVals, 10^-2.5]]; 
segs =  BreakAtXVal[Apply[Join, ToFunctions /@ 
  (elimDups2DOrdered[#, 10^-4] & /@ jj)], criticalXVals];
 Table[vv = Select[Select[segs, 
               Length[elimDups1D[ First /@ # , 10^-2.5]] != 1 & ], 
 (* eliminates near-vertical segments *)
   Abs[#[[1, 1]] - criticalXVals[[i]]] <= 10^-2.5 &];
       xxx = (First @ Thread[#] &) /@ vv;
       mmm = Append[MapIndexed[
                (vv[[First[#2], #1]] + vv[[First[#2], #1+1]] )/2 & , 
                   (home1[#, 
                    (criticalXVals[[i]] + criticalXVals[[i+1]])/2] & /@ xxx)],
 
               {(criticalXVals[[i]] + criticalXVals[[i+1]])/2, N@ymax}];    
      
      (* here we are matching the segments with their evaluation value, 
      for ease and correctness of sorting *)
      vvv = Drop[First /@ Sort[Thread[
           {Append[vv, {}], mmm}], #1[[-1, -1]]  < #2[[-1, -1]]  &], -1];
      mmm = Sort[mmm, #1[[-1]] < #2[[-1]] &];
      grays = GrayLevel /@ (( (Sign[N@funcs] /. Thread[{x,y}->N[#]]) /. {
        {1,1}->.32, {0, 0} -> 1, {1, -1}->.6, {-1,-1}->.46, {-1,1}->.8,
        {0, _} :> 1, {_, 0} :> 1 })
        & /@   ((Drop[mmm, -1] + Rest[mmm])/2)   );
      
Thread[{grays, Polygon[Join[#, {{#[[-1,1]], ymax}, 
           {#[[1,1]], ymax}}]]& /@ vvv}],
  {i, 1, Length[criticalXVals]-1}]}]], {} ],   

Graphics[{If[showLines, Join[cs /. Automatic->{},
   {Append[nullcSty[[1]], Line /@ contourData[[1]]]},
   {Append[nullcSty[[2]], Line /@ contourData[[2]]]}], {}],

If[stabQ, Append[stabCol, Point /@ stabPts], {}]}], 
FilterOptions[Graphics, opts]]]


(* ****************** PhaseLine  *************  *)

Options[PhaseLine] = {
  FieldColor->GrayLevel[0],
  NumberFish->10,
  FieldThickness->Thickness[0.004],
  AspectRatio->0.1,
  TimeScale->10,
  FlowField->False,
  FieldLogScale->1,
  FlowThickness->1};

PhaseLine[e_?FreeOfSetQ, x___] :=  (
    Message[VDS::insarg, PhaseLine, 3]; 
    HoldForm[PhaseLine[e,x]]) /; Length[{x}] < 2

PhaseLine[e_?FreeOfSetQ, x___] :=  (Message[VDS::nonopt, 3]; 
      HoldForm[PhaseLine[e, x]]) /; 
     Length[{x}] != 2 && Union[Head /@ Drop[{x}, 2]] =!= {Rule}

PhaseLine[e_?FreeOfSetQ, x___] :=  (Message[VDS::badarg]; 
            HoldForm[PhaseLine[e,x]]) /; Head[e] =!= Equal

PhaseLine[e_?FreeOfSetQ, x___] :=  (Message[VDS::badarg, 3, List]; 
            HoldForm[PhaseLine[e,x]]) /; !ListQ[{x}[[2]]]

PhaseLine[e_?FreeOfSetQ, x___] :=  (Message[VDS::nonnum, {3,2}]; 
           HoldForm[PhaseLine[e,x]]) /; !NumberQ[N[{x}[[2,2]]]]

PhaseLine[e_?FreeOfSetQ, x___] :=  (Message[VDS::nonnum, {3,3}]; 
           HoldForm[PhaseLine[e,x]]) /; !NumberQ[N[{x}[[2,3]]]]

PhaseLine[e_?FreeOfSetQ, x___] :=  (Message[VDS::badarg, {3,1}]; 
           HoldForm[PhaseLine[e,x]]) /; !({x}[[2, 1]] === Head[{x}[[1]]]) 


PhaseLine[equation_?FreeOfSetQ, unk_, {x_, tx0_, tx1_}, opts___] :=
 Module[{t, t0, t1}, 
 t = unk[[-1]];
{col, fhead, thick, fish, asp, plotR,
     numfish, StartPts, ttt, fishl} = 
{FieldColor, FlowThickness, FieldThickness, FlowField, AspectRatio, 
     PlotRange, NumberFish, Steps, TimeScale, FieldLogScale} /. {opts} /. 
     Options[PhaseLine] /. Options[FindAllCrossings];

{x0, x1} = N[{tx0, tx1}]; 
testx = (x1 - x0)/Max[{20, StartPts}];          
xinterval = Interval[{x0,x1}];

equil = Append[Prepend[FindAllCrossings[equation[[2]] /. x[t] -> x, 
       {x, x0, x1}], x0], x1];       
       
grayareas = Map[{GrayLevel[(ccol = 
    Sign[N[equation[[2]] /. x[t] -> ((bb = #[[1]]) +
         Random[]*((cc = #[[2]])-bb))]]; 
         If[ccol === 0, 1, (ccol + 4)/6] )],
    Polygon[{{bb, -(aa = 1)}, {bb, 1}, {cc, 1}, {cc, -aa}}]}&,
         Drop[Transpose[{equil, RotateLeft[equil]}], -1]];

{bx, by} = {(x1 - x0) asp, 3}/40*fhead;
dt = ttt/numfish;

makefish[{si_, sf_}] := If[fish,
   {col, Polygon[{{si, 0}, {sf, 2 by}, {sf, -2 by}}],
        Disk[{sf, 0}, 2*{bx, by}, {0, 2Pi}]},      
   {thick, col, Line[{{si, 0}, {sf, 0}, {sf - 
        (bbx = 8 If[sf == si, 0, Sign[sf - si]*Exp[Log[fishl/10]*
         Log[10, maxfish/(Abs[sf - si])]]])*
    bx, bbx*by}, {sf, 0}, {sf - bbx*bx, -bbx*by}}]}];      
      
allpoints = {};
StartPts = Range[x0, x1, (x1 - x0)/StartPts];
xdone = Interval[{x0 - 1, x0 - 0.5}];  
idex = -1;
   
While[StartPts =!= {}, (  
idex *= -1;
xx0 = N[StartPts[[idex]]]; 
If[(xini = N[equation[[2]] /. x[t] -> xx0]) == 0. || !NumberQ[xini],
    (allpoints =   Join[allpoints, {{xx0,xx0}}];
       StartPts = Drop[StartPts,idex]), (

If[idex*xini > 0,
    {tt0, tt1} = N[{0, ttt}],
    {tt0, tt1} = N[{-ttt, 0}]];

sol = First[x[t] /. NDSolve[
  {equation, x[0] == xx0}, x[t], 
     {t, tt0, tt1}, Evaluate[FilterOptions[NDSolve, opts]]]];
solf = sol[[0]]; 
{t0, t1} = getDomain[sol]; 

If[t1 > t0, (tt1 = 0;
  If[t1 > 0, 
  Scan[If[(IntervalMemberQ[xinterval, (xaaa = #[[2]])]) &&
         !IntervalMemberQ[xdone, xaaa],
         tt1 = #[[1]],
         Null]&, GetPtsFull[sol]], 
  Scan[If[(IntervalMemberQ[xinterval, (xaaa = #[[2]])]) &&
         !IntervalMemberQ[xdone, xaaa],
         tt1 = #[[1]],
         Null]&, Reverse[GetPtsFull[sol]] ]];
         
xdone = IntervalUnion[xdone, 
   Interval[{sss0 = Re[solf[t0]], sss1 = Re[solf[t1]]}],
   Interval[{sss0 - testx/2, sss0 + testx/2}],
   Interval[{sss1 - testx/2, sss1 + testx/2}]];
   
If[t1 > 0, t1 = tt1, t0 = tt1];
StartPts = DeleteCases[Map[If[IntervalMemberQ[
    xdone, #] , Null, #]&, StartPts], Null];
    
points = Map[{solf[#], solf[# + 4dt/5]}&, 
            Range[t0,  t1 - 4dt/5, dt]]; 
           
allpoints = Join[allpoints, points]),
  StartPts = Drop[StartPts,idex]])])];

maxfish = Max[Map[Abs[#[[1]] - #[[2]]]&, allpoints]];
fish = Map[makefish, allpoints];

Show[Graphics[{grayareas, fish}],
  Evaluate[FilterOptions[Graphics, opts]],   
  Frame -> True, Axes->None,
  FrameTicks -> {Automatic, None, None, None},
  FrameStyle -> {{GrayLevel[0]}, {GrayLevel[1]}, {GrayLevel[1]}, {GrayLevel[
1]}},
  AspectRatio -> asp, 
 (*  Axes -> {True, False}, *)
  PlotRange -> {{x0, x1}, {-2, 1}}
(*,  AxesOrigin -> {x0, -1.2} *) ]]


(*  *****************  PhasePlot **************** *)
 
 Options[PhasePlot] = {
     AspectRatio->1/GoldenRatio,
     Background->GrayLevel[1],
     BoxRatios->{1, 1, 0.4},
     ComputeWindow -> Automatic,
     DefaultColor->Automatic,
     DirectionArrow->False,
     EquilibriumPointStyle->{GrayLevel[0], PointSize[0.021]},
     FastPlotting->False,
     FieldMeshSize->15,
     FlowField->False,
     Frame->False,
     InitialPointStyle->{GrayLevel[0], PointSize[0.015]},
     InitialValues->{},
     Method->Automatic,
     Nullclines->False,
     NullclineShading->False,
     PlotStyle->{},
     ParametricPlotFunction->ParametricPlot,
     PrecisionGoal->Automatic,
     Rainbow->False,
     SaveLastPoint->False,
     SaveSolution->False,
     ShowInitialValues->False,
     ShowNumberOfSteps -> False,
     ShowOrbits->True,
     ShowEquilibria->False,
     SolutionName->"Solution",
     Source->FirstOrder,
     StayInWindow->False,
     Ticks->Automatic,
     TickLabelsOnly->False,
     VectorField->False,
     WindowShade->None,
     WorkingPrecision->$MachinePrecision};     
     
(* case of no t-iterator *)
PhasePlot[eqn_?FreeOfSetQ /; Length[eqn] == 2, unks_,
  {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, opts___Rule] := If[
    (InitialValues /. {opts} /. Options[PhasePlot]) =!= {},
     Message[PhasePlot::incons]; HoldForm[PhasePlot[eqn, unks,
       {x, xmin, xmax}, {y, ymin, ymax}, opts]], 
  PhasePlot[eqn, unks, {unks[[1, -1]], 0, 1}, {x, xmin, xmax}, 
  {y, ymin, ymax}, opts]]
  
PhasePlot[e_?FreeOfSetQ, x___] :=  (Message[VDS::insarg, PhasePlot, 4]; 
    HoldForm[PhasePlot[e,x]]) /; Length[{x}] < 3

PhasePlot[e_?FreeOfSetQ, x__] :=  (Message[VDS::badarg]; 
            HoldForm[PhasePlot[e,x]]) /; 
  (Position[Head /@ Take[{x}, 3], _?(# =!= List &)]) != {{}}

PhasePlot[e_?FreeOfSetQ, x__] :=  (Message[VDS::badarg]; 
            HoldForm[PhasePlot[e,x]]) /; NumberQ[N[{x}[[1,2]]]]

PhasePlot[e_?FreeOfSetQ, x__] :=  (Message[VDS::nonnum, {r[[1,1]]+3, 2}]; 
            HoldForm[PhasePlot[e,x]]) /;
        (r=Position[{x}[[{2,3}, 2]], 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}

PhasePlot[e_?FreeOfSetQ, x__] :=  (Message[VDS::nonnum, {r[[1,1]]+3, 3}]; 
    HoldForm[PhasePlot[e,x]]) /;  (r=Position[{x}[[{2,3}, 3]], 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}

(* case of t-iterator and 2 display variables *)

PhasePlot[equation_?FreeOfSetQ /; Length[equation] > 1, unks_,
  {t_, tmin_, tmax_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, opts___Rule] :=
Module[{},
num = Length[unks];
halt = HoldForm[PhasePlot[equation, unks, {t, tmin, tmax}, 
                  {x, xmin, xmax}, {y, ymin, ymax}, opts]];

If[!ListQ[equation], Message[VDS::notlist]; Return[halt]];

If[!CheckOptions[{PhasePlot, Graphics, NDSolve, RungeKutta4, Euler, 
  ContourPlot, Plot, DirectionFieldPlot, Equilibria, ColorParametricPlot,
  ProjectionPlot3D, Graphics3D, NullclinePlot, FlowParametricPlot}, opts], 
    Return[halt]];      

If[Or @@ (FreeQ[equation, Head[#]'[t]] & /@ unks), Message[VDS::badform, t]; 
  Return[halt]];

(* the following line sets newOpts to be all the options that should
be used, whether they come from Defaults or were added in via opts.
It also fixes a bug in DefaultColor *)

newOpts = Sequence @@ (Join[
  Options[PhasePlot] /. ((#[[1]]->_) :> # &) /@ {opts},
  Complement[{opts}, {FilterOptions[PhasePlot, opts]}]] /. 
       RGBColor[a_, a_, a_] :> GrayLevel[a]);

{fieldSteps, ptsQ, fieldQ, orbSty, pts, ptCol, backCol, asp, ppfcn, saveSol,
 lastQ, arr, stableQ, stabCol, nullclinesQ, grayregQ, showOrbsQ, met, 
 showFish, stayIn, defC, fastPlQ, sou, solName, tic, ticLO, compWin, winShad,
 rain, wp, showNumQ} =   {FieldMeshSize, ShowInitialValues, VectorField, 
PlotStyle, 
 InitialValues, InitialPointStyle, Background, AspectRatio,
 ParametricPlotFunction, SaveSolution, SaveLastPoint, 
 DirectionArrow, ShowEquilibria, EquilibriumPointStyle, Nullclines, 
 NullclineShading, ShowOrbits, Method, FlowField, StayInWindow,
 DefaultColor, FastPlotting, Source, SolutionName, Ticks,
 TickLabelsOnly, ComputeWindow, WindowShade, Rainbow,
 WorkingPrecision, ShowNumberOfSteps} /. {newOpts};     

eqn = If[sou =!= SecondOrder, NormalForm[equation, unks], equation];

method = met /. Automatic->NDSolve;
saveSol = saveSol || MemberQ[First /@ {opts}, SolutionName];
If[!ListQ[stabCol], stabCol = {stabCol}];
If[!ListQ[ptCol], ptCol = {ptCol}];

xminn = Min[xmin, xmax]; xmaxx = Max[xmin, xmax];
ymaxx = Max[ymin, ymax]; yminn = Min[ymin, ymax];

{{cwx0, cwx1}, {cwy0, cwy1}} = 
  N[compWin /. Automatic -> {{xminn, xmaxx}, {yminn, ymaxx}}];
    
Off[NDSolve::ndnum];

If[method =!= NDSolve && stayIn, Message[VisualDSolve::bdmet]; Return[halt]];

c1 = (First /@ eqn) /. xx_'[t] :> xx;
c2 = {x, y};
c3 = Head /@ unks;
If[(Select[c1, MemberQ[c2, #]&] =!= c2)  || c1 =!= c3,
   Message[PhasePlot::bdordr]; Return[halt]];
   
If[!StringQ[solName], Message[VisualDSolve::badname]; Return[halt]];
rhSides = (Last /@ eqn) /. {x[t]->x, y[t]->y};
axLab = AxesLabel /. {newOpts} /. AxesLabel->None;
arrows = {};
If[saveSol, solNameTemp = {}];
If[lastQ, LastPoint = {}];
If[stayIn, WindowExitData = {}];
projQ = (ppfcn === ProjectionPlot3D);
font[s_] := FontForm[s, {"Courier", size}];
delY = 2 10/72 (ymaxx - yminn) / 5;
pRange = N[{{xminn, xmaxx}, {yminn, ymaxx + If[axLab =!= None, 2 delY, 0]}}];
    
orbSty = fixSty[orbSty, Length[pts]];
{xR, yR} = {xmaxx - xminn, ymaxx - yminn};

If[pts != {} && NumberQ[N[pts[[1]]]], pts = {pts}];
If[pts == {} && !fieldQ && !showFish && !nullclinesQ && !grayregQ && !stableQ,

  Message[PhasePlot::nothng]; Return[halt]];
defC = defC /. Automatic->(backCol /. z_?NumberQ :> 1. - z);
asp = N[asp /. Automatic->yR/xR];

If[NumberQ[fieldSteps], fieldSteps = {fieldSteps, fieldSteps}];

If[pts != {}, pts = N[Map[If[NumberQ[N[#[[-1]]]], 
  Append[#, {tmin, tmax}], #]&, pts], If[wp > 16, wp, 15]]];
If[pts != {}, pts = N[Map[
If[Length[#] != num+2, Prepend[#, tmin], #]&, pts], If[wp > 16, wp, 15]]]; 
transformed = False;
nSoln[ut_, xVals__, {t1_, t2_}] := (
 
 If[stayIn && N[{t1, t2}] == N[{tmin, tmax}] && !transformed &&
     !ver3Q, 
  transformed = True;
  eqn = ReplacePart[eqn, eqn[[1,1]] == (If[
       cwx0 <= x[t] <=  cwx1 && 
        cwy0 <= y[t] <= cwy1,
          Evaluate[eqn[[1, 2]]]]), 1]];

stopOpt = If[stayIn && ver3Q, 
    StoppingTest -> !(cwx0 <= x[t] <=  cwx1 && 
        cwy0 <= y[t] <= cwy1), WorkingPrecision -> wp];

First[If[Length[pts] > 1,
Check[stemp = method[
    Join[eqn, (#[[1]] /. t->ut) == #[[2]] & /@ Transpose[{unks, {xVals}}]], 
    #[t] & /@ (First /@ Head /@ First /@ eqn),
    {t, If[stayIn && {t1, t2} == {tmin, tmax}, 2 t1 - t2, t1], t2}, 
    FilterOptions[method, newOpts], Evaluate[stopOpt]],
  Message[VDS::mxstep, {ut, xVals}, t,  getDomain[stemp[[1,1,2]]][[2]]]; 
stemp,
 NDSolve::mxst],
 
 method[Join[eqn, (#[[1]] /. t->ut) == #[[2]] & /@ Transpose[{unks, {
xVals}}]], 
   #[t] & /@ (First /@ Head /@ First /@ eqn),
   {t, If[stayIn && {t1, t2} == {tmin, tmax}, 2 t1 - t2, t1], t2}, 
   FilterOptions[method, newOpts], Evaluate[stopOpt]]]]);

 If[!FreeQ[rhSides, t] && pts == {}, 
     Message[PhasePlot::noinv]; Return[halt]];

If[(fieldQ || showFish) && !FreeQ[rhSides, t], 
  Message[PhasePlot::ntatn]; Return[halt]];

Catch[If[stayIn, {#, WindowExitData}, #]& @
Show[

Graphics[If[winShad === None, {}, 
{{winShad, Rectangle[{xminn , yminn }, {xmaxx, ymaxx}]}}]],
   
If[(nullclinesQ || grayregQ) && num == 2 && FreeQ[rhSides, t], 
   NullclinePlot[rhSides,  
     {x, xmin, xmax}, {y, ymin, ymax},
      DisplayFunction->Identity,
      Nullclines->nullclinesQ,
      ShowEquilibria->stableQ,
      FilterOptions[NullclinePlot, opts],
      FilterOptions[Equilibria, opts]], {}],

If[fieldQ || showFish, 
   N[DirectionFieldPlot[rhSides,
   {x, xminn, xmaxx, xR/fieldSteps[[1]]}, 
   {y, yminn, ymaxx, yR/fieldSteps[[2]]},
   Axes->None,AxesLabel->None, Source->PhasePlot,
   Evaluate[Sequence[FilterOptions[DirectionFieldPlot, newOpts], 
                     FilterOptions[Graphics, newOpts]]],
   PlotRange->pRange,
   DisplayFunction->Identity]], {}], 

If[!showOrbsQ || pts == {}, {}, MapIndexed[
   (sol = nSoln @@ #1;
 If[method === NDSolve && showNumQ, Print[StringForm[
   "The initial value `` required `` steps.", 
      Drop[#1, -1], Length[GetPts[sol[[1,2]]]]]]];        
iter = If[N[#1[[-1]]] != N[{tmin, tmax}],
           {t, Max[getDomain[sol[[1,2]]][[1]], #1[[-1, 1]]], 
               Min[getDomain[sol[[1,2]]][[2]], #1[[-1, 2]]]},
           {t, getDomain[sol[[1,2]]][[1]], getDomain[sol[[1,2]]][[2]]}];  


If[arr,
      {temp1, temp2} =  {x[t], y[t]} /. sol /. t->(iter[[2]]+iter[[3]])/2;
      {slope1, slope2} =
           (Last /@ (eqn[[Flatten[ Position[unks, #] & /@ {x[t], y[t]}]]])) /.

          {x[t]->temp1, y[t]->temp2, t->(iter[[2]]+iter[[3]])/2};

{bx, by} = If[asp<1, {(pRange[[1,2]]-pRange[[1,1]]),
                      (pRange[[2,2]]-pRange[[2,1]])/asp},
             {asp*(pRange[[1,2]]-pRange[[1,1]]),
              (pRange[[2,2]] - pRange[[2,1]])}]/25;
             
arrowtip = {slope1, slope2} * (gg = g[bx, by, slope1, slope2]);
arrowedge = {-slope2 * bx / by, slope1 * by / bx} *gg / 2.5;
AppendTo[arrows, {defC, Polygon[{temp1, temp2}  + # & /@      
  {0, -arrowtip/4 + arrowedge, arrowtip, -arrowtip/4 - arrowedge}   ]}]];
          
  If[saveSol, AppendTo[solNameTemp, {x[t], y[t]} /. sol]];
  
  If[stayIn, AppendTo[WindowExitData, 
    ReplacePart[#1,  getDomain[sol[[1,2]]], -1]]];
  If[lastQ, AppendTo[LastPoint, unks /. sol /. t -> 
  getDomain[sol[[1,2]]][[2]]  ]];
  
  If[!(And @@ ( {x[t], y[t]} /. sol /. 
      InterpolatingFunction[___][t] :> True)),
       Message[VisualDSolve::failure]; Graphics[{}],
 
If[fastPlQ && ppfcn === ParametricPlot,
    Graphics[Append[If[rain,  
       Prepend[orbSty[[First[#2]]], Hue[1 - 0.8 #2[[1]] / Length[pts]]],
        orbSty[[First[#2]]]],   
    dompts = getDomPts[x[t] /. sol];
    ic = 1;
    While[ (!stayIn && dompts[[ic]] < Min[#1[[-1]]]) ||
           (stayIn && dompts[[ic]] < 2 Min[#1[[-1]]] - Max[#1[[-1]]] ), ic++];

    jc = ic;
While[jc <= Length[dompts] && dompts[[jc]] <= Max[#1[[-1]]], jc++];
Line[Take[Thread[GetPts /@ ({x[t], y[t]} /. sol)], {ic, jc-1}]]]],

 If[iter[[2]] < iter[[3]], 
 
  If[ppfcn === ProjectionPlot3D, Throw[Return[#]]&, Identity][
   pptemp = 
   If[ppfcn === FlowParametricPlot, 
     FlowParametricPlot[{x[t], y[t]} /. sol, {t, tmin, tmax},
       {x, xminn, xmaxx}, {y, yminn, ymaxx},
       Evaluate[FilterOptions[FlowParametricPlot, 
        Sequence @@ ({newOpts} /. ((Rainbow -> True) -> 
         (FlowColor -> Hue[1-0.8*#2[[1]]/Length[pts]])
            ))]]], 
   
   ppfcn[Evaluate[{x[t], y[t]} /. sol], Evaluate[iter],   
   Evaluate[Sequence @@ (If[projQ, {{xminn, xmaxx}, {yminn, ymaxx}}, {}])],
   PlotStyle->Evaluate[If[rain,
       {Prepend[orbSty[[First[#2]]], Hue[1 - 0.8 #2[[1]] / Length[pts]]]},
        orbSty[[#2]]]],
   Evaluate[FilterOptions[ppfcn, newOpts]],
   Evaluate[FilterOptions[If[projQ, Graphics3D, ppfcn], newOpts]],
      DisplayFunction->If[projQ, $DisplayFunction, Identity]]];
    If[Head[pptemp] === HoldForm, 
   Throw[Return[#]]&, Identity][pptemp]], {}]] ])&, pts]],
  
On[NDSolve::ndnum];
solNameTemp = DeleteCases[DeleteCases[solNameTemp, 
  InterpolatingFunction[{a_, a_} | {{a_, a_}}, b__][t], 2], {}];
If[Length[solNameTemp] == 1, solNameTemp = First[solNameTemp]];

Clear[Evaluate@StringJoin["Global`", solName]];
Evaluate[ToExpression@StringJoin["Global`", solName]] = solNameTemp;
Clear[solNameTemp];

If[Length[LastPoint] == 1, LastPoint = First[LastPoint]];  

Graphics[{
  If[ptsQ && pts != {}, Join[ptCol,  
    Point[Rest[Drop[#, -1]][[
       Flatten[Join[Position[unks, x[t]], Position[unks, y[t]]]]]]] & /@ 
       pts], {}],
  
  
    If[stableQ && (Length[eqn] == 2), Append[stabCol, 
     Point /@ If[nullclinesQ || grayregQ, eqPts,
        Equilibria[Last /@ (eqn /. {x[t]->x, y[t]->y}), 
          {x, xmin, xmax}, {y, ymin, ymax}, 
          Evaluate[FilterOptions[Equilibria, opts]],
          
          EquilibriumFindingMethod -> If[And @@ (FreeQ[eqn, #] & /@ {Abs, 
Sign, If, Which}),
            FindRoot, FindMinMethod] ]]],   {}],
 arrows,

 If[axLab =!= None,{
     {backCol, Rectangle[{xminn - xR/100, ymaxx}, 
        {xmaxx + xR/100, pRange[[2,2]]}]},
      Text[axLab[[2]], {xminn, ymaxx + yR/50}, {-1, -1}]}, {}],
 
{Thickness[0.004], defC,
  Line[{{xmin, ymin}, {xmin, ymax}, {xmax, ymax}, {xmax, ymin},
          {xmin,ymin}}]}}],
 DefaultColor->defC,
 AxesLabel->If[(axLab =!= None),
   {axLab[[1]], None}, axLab],
 Sequence @@ DeleteCases[{FilterOptions[Graphics, newOpts]},
   Ticks -> _],
 Ticks -> (If[grayregQ, fixTicks, Identity][
    If[tic === Automatic, (
    Map[Flatten, 
    FullOptions[Graphics[{Line[{{xminn, yminn}, {xmaxx, ymaxx}}]},
     Axes->Automatic, AxesOrigin->{xminn, yminn}, Ticks -> Automatic], Ticks]
 /. 
        {a_, b_, _, _} :> If[b == "", {}, a]]), tic]]), 
 Axes->Automatic,
 AxesOrigin->{xmin, ymin},
 AxesStyle->{backCol, AbsoluteThickness[1]},
 PlotRange->pRange,
 DisplayFunction->Evaluate[If[projQ, Identity, $DisplayFunction]]]] ]



 (*  *****************  PhasePlot, 3D case **************** *) 

(* this option list is for internal use only! *)

Options[PhasePlot3DVersion] = {
     ShowInitialValues->False,
     PlotPoints->100,
     FastPlotting -> False,
     DefaultColor->Automatic,
     PrecisionGoal->Automatic,
     PlotStyle -> {AbsoluteThickness[0.5]},
     AspectRatio->1/GoldenRatio,
     InitialValues->Automatic,
     Axes->True,
     InitialPointStyle->{GrayLevel[0], PointSize[0.021]},
     Ticks->Automatic,
     Background->GrayLevel[1],
     BoxRatios->{1, 1, 0.4},
     ParametricPlotFunction->ParametricPlot3D,
     Rainbow->False,
     SaveSolution->False,
     ShowNumberOfSteps -> False,
     SolutionName->"Solution",
     SaveLastPoint->False,
     WorkingPrecision->$MachinePrecision};
 
Off[ParametricPlot::ppcom];

PhasePlot[eqn_?FreeOfSetQ, unks_List,
 {t_, tmin_, tmax_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_},
 {z_, zmin_, zmax_}, opts___Rule]  :=
  
 Module[{nSoln},
halt = HoldForm[PhasePlot[eqn, unks, {t, tmin, tmax}, 
    {x, xmin, xmax}, {y, ymin, ymax}, {z, zmin, zmax}, opts]];

If[!ListQ[eqn], Message[VDS::notlist]; Return[halt]];

If[!CheckOptions[{PhasePlot3DVersion, NDSolve, ParametricPlot3D}, opts],
  Return[halt]];
      
(* Because SetOptions will not affect this 3D version, we fix now,
for the most common type of SetOption use only *)

If[FastPlotting /. Options[PhasePlot], 
  PrependTo[Options[PhasePlot3DVersion], FastPlotting -> True]];
            
(* the following line sets newOpts to be all the options that should
be used, whether they come from Defaults or were added in via opts *)

newOpts = Apply[Sequence,
 Join[{PlotPoints->(PlotPoints /. {opts} /. PlotPoints->100)},
 Options[PhasePlot3DVersion] /. ((#[[1]]->_) :> # &) /@ {opts},
 Complement[{opts}, {FilterOptions[PhasePlot3DVersion, opts]}] ]];

{showNumQ, numberPts, ptsQ, orbSty, pts, ptCol, backCol, asp, 
 saveSol, lastQ, ppFcn, defC, solName, rain, wp, fastQ} = 
 {ShowNumberOfSteps, PlotPoints, ShowInitialValues, PlotStyle, InitialValues, 

 InitialPointStyle, Background, AspectRatio, SaveSolution,
 SaveLastPoint, ParametricPlotFunction,
 DefaultColor, SolutionName, Rainbow, WorkingPrecision,
 FastPlotting} /. {newOpts};
 
eqn1 = NormalForm[eqn, unks];  (* to fix 2.2 bug *)

nSoln[{ut_, xVals__, {t1_, t2_}}] := First[NDSolve[
  Join[eqn1, (#[[1]] /. t->ut) == #[[2]] & /@ Transpose[{unks, {xVals}}]], 
  First /@ Head /@ First /@ eqn1, {t, t1, t2}, 
  FilterOptions[NDSolve, newOpts]]];
   
saveSol = saveSol || MemberQ[First /@ {opts}, SolutionName];
c1 = Cases[eqn1, xx_'[_], Infinity] /. zz_'[t] :> zz;
c2 = {x, y, z};
c3 = Head /@ unks;
If[(Select[c1, MemberQ[c2, #]&] =!= c2)  || c1 =!= c3,
    Message[PhasePlot::bdordr]; Return[halt]];

If[!StringQ[solName], Message[VisualDSolve::badname]; Return[halt]];
solTemp = {};
LastPoint = {};
{xR, yR, zR} = {xmax - xmin, ymax - ymin, zmax - zmin};
If[TensorRank[pts] == 1, pts = {pts}];
defC = defC /. Automatic->(backCol /. zz_?NumberQ :> 1. - zz);
asp = asp /. Automatic->yR/xR;
axLab = AxesLabel /. {newOpts} /. AxesLabel->None;

pts = N[pts /. 
  {u__?(NumberQ[N[#]]&), {t1_?(NumberQ[N[#]]&), t2_?(NumberQ[N[#]]&)}} :> 
                 {u, {t1, t2, xx}} /.
     {u__?(NumberQ[N[#]]&)} :> {u, {tmin, tmax}} /.
     {u__, {t1_, t2_, xx}} :> {u, {t1, t2}} /.
  (Automatic->{Append[(unks /.
           {x[t]->(xmin+xmax)/2, y[t]->(ymin+ymax)/2, z[t]->(zmin+zmax)/2} /.
   (#->0 & /@ Complement[unks, {x[t], y[t], z[t]}])), {tmin, tmax}]}), wp];
   
If[pts != {}, pts = N[Map[
  If[Length[#] != Length[unks]+2, Prepend[#, tmin], #]&, pts], If[wp > 16, wp,
 15]]];
   
If[TensorRank[pts] == 3, pts = pts[[1]]];

orbSty = fixSty[orbSty, Length[pts]];
If[!ListQ[ptCol], ptCol = {ptCol}];
 
Show[
   MapIndexed[
   (sol = nSoln[#1]; 
    If[showNumQ, Print[StringForm[
    "The initial value `` required `` steps.", 
       Drop[#1, -1], Length[GetPts[sol[[1, 2]][t]]]]]];
   
    If[saveSol, AppendTo[solTemp, {x[t], y[t], z[t]} /. sol]];
    If[lastQ, AppendTo[LastPoint, unks /. sol /. t->#1[[-1, 2]]]];
    If[!(And @@ ( {x[t], y[t], z[t]} /. sol  /. 
       InterpolatingFunction[___][t] :> True)),
       Message[VisualDSolve::failure]; Graphics[{}],
  
  fffx = Compile[{t}, Evaluate[x[t] /. sol]];
  fffy = Compile[{t}, Evaluate[y[t] /. sol]]; 
  fffz = Compile[{t}, Evaluate[z[t] /. sol]]; 
 Show[If[ppFcn === ParametricPlot3D,
   Graphics3D[Append[
     If[rain,  
       Prepend[orbSty[[First[#2]]], Hue[1 - 0.8 #2[[1]] / Length[pts]]],
        orbSty[[First[#2]]]],
   If[fastQ, Line[Thread[Map[Function[ll, GetPts[ll, #1[[-1]]]], {x[t], y[t], z[t]} /. sol]]],  
    Line[Table[{fffx[Min[tt, tmax]], fffy[Min[tt, tmax]], fffz[Min[tt, 
tmax]]},         
      {tt, #1[[-1, 1]], #1[[-1, 2]], 
         (#1[[-1, 2]] - #1[[-1, 1]])/(numberPts-1)} ]]]],
     Evaluate[FilterOptions[Graphics3D, newOpts]], 
       DisplayFunction->Identity],
   
    (* else *) ColorParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], 
     {t, #1[[-1, 1]], #1[[-1, 2]]},
    PlotStyle->orbSty, DisplayFunction->Identity,
   Evaluate[FilterOptions[ColorParametricPlot3D, newOpts]] ]]]])&, pts],

If[Length[solTemp] == 1, solTemp = First[solTemp]];
If[Length[LastPoint] == 1, LastPoint = First[LastPoint]];

Clear[Evaluate@StringJoin["Global`", solName]];
Evaluate[ToExpression@StringJoin["Global`", solName]] = solTemp;
Clear[solTemp];

Graphics3D[{
  If[ptsQ, Append[ptCol, 
     Point[#[[Flatten[
       Join[Position[unks, x[t]], Position[unks, y[t]], 
       Position[unks, z[t]]]]]]] & /@ 
       pts], {}]   }],

 Sequence @@
  DeleteCases[ {FilterOptions[Graphics3D, newOpts]}, AxesLabel->_],
Axes->Automatic, Boxed->False,
AxesLabel-> (ToString /@ {x, y, z}), 
DefaultColor->defC,
PlotRange->{{xmin, xmax}, {ymin, ymax}, {zmin, zmax}},
 DisplayFunction->$DisplayFunction]]




(* ************ PoincareSection ************ *)

Options[PoincareSection] = {TimeInterval->Automatic, 
    WorkingPrecision->$MachinePrecision,
    AccuracyGoal->Automatic,
    PrecisionGoal->Automatic,
    PointStyle->{PointSize[0.004]},
    Rainbow -> False,
    SavePoints -> False,
    WindowShade->None};

PoincareSection[x___] :=  (Message[VDS::insarg, PoincareSection, 5]; 
    HoldForm[PoincareSection[x]]) /; Length[{x}] < 5

PoincareSection[x__] :=  (
  Message[VDS::nonopt, 5]; HoldForm[PoincareSection[x]]) /; 
     Length[{x}] != 5 && Union[Head /@ Drop[{x}, 5]] =!= {Rule}

PoincareSection[x__] :=  (Message[VDS::badarg]; 
            HoldForm[PoincareSection[x]]) /; 
        (Position[Head /@ Rest[Take[List[x], 5]], _?(# =!= List &)]) != {{}}

PoincareSection[x__] :=  (Message[VDS::nonnum, {r[[1,1]]+1, 2}]; 
            HoldForm[PoincareSection[x]]) /;
        (r=Position[(List[x][[{2,3,4}, 2]]), 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}
PoincareSection[x__] :=  (Message[VDS::nonnum, {r[[1,1]]+1, 3}]; 
            HoldForm[PoincareSection[x]]) /; 
        (r=Position[(List[x][[{2,3,4}, 3]]), 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}
PoincareSection[x__] :=  (Message[VDS::nonnum, {5, r[[1,1]]}]; 
            HoldForm[PoincareSection[x]]) /; 
        (r=Position[(List[x][[5, {1,2}]]), 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}

PoincareSection[{eq1_, eq2_}, {t_, tmin_, tmax_}, {x_, xmin_, xmax_}, 
  {y_, ymin_, ymax_}, {x0_, y0_}, opts___] := Module[{ptss},
   {delt, sty, ag, pg, wp, winShad, rnbow, showPts} = {TimeInterval, 
PointStyle,
  AccuracyGoal, PrecisionGoal, WorkingPrecision,
  WindowShade, Rainbow, SavePoints} /. {opts} /. Options[PoincareSection];
       
halt = HoldForm[PoincareSection[{eq1, eq2}, {t, tmin, tmax}, {x, xmin, xmax}, 

  {y, ymin, ymax}, {x0, y0}, opts]];
If[!NumberQ[N@delt] && delt =!= Automatic, 
      Message[VDS::badopt, delt, TimeInterval]; Return[halt]];
If[!IntegerQ[wp] && wp =!= Automatic,
   Message[VDS::badopt, wp, WorkingPrecision]; Return[halt]];

If[!IntegerQ[ag] && ag =!= Automatic, 
   Message[VDS::badopt, ag, AccuracyGoal]; Return[halt]];

If[!IntegerQ[pg] && pg =!= Automatic, 
    Message[VDS::badopt, pg, PrecisionGoal]; Return[halt]];

If[!MemberQ[{RGBColor, GrayLevel, Hue}, Head[winShad]] && winShad =!= None, 
    Message[VDS::badopt, winShad, WindowShade]; Return[halt]];
 
delt = delt /. Automatic->(tmax-tmin)/100;
wp = wp /. Automatic->$MachinePrecision;
If[!ListQ[sty] && sty =!= Rainbow, sty = {sty}];

mxChk = Check[NDSolve[
  Join[{eq1, eq2}, {x[tmin] == x0, y[tmin] == y0}], 
         {x[t], y[t]}, {t, tmin, tmax}, FilterOptions[NDSolve, opts], 
              MaxSteps -> 1000],
   Message[PoincareSection::mxst]; False, NDSolve::mxst];
If[mxChk === False, Return[halt]];

pts1 = Thread[
First[{x[t], y[t]} /. mxChk] /. t->N[Range[tmin, tmax, delt], wp]];
ptss = Point /@ pts1;
If[rnbow, ptss = 
  Transpose[{Hue /@ Range[0, 0.5, 0.5/(Length[ptss]-1)], ptss}]];
  
Show[gg=Graphics[{

If[winShad === None, {}, 
  {{winShad, Rectangle[{xmin, ymin}, {xmax, ymax}]},
    Thickness[0.004], Line[{{xmin, ymin},
       {xmin, ymax},{xmax, ymax},{xmax, ymin}, {xmin, ymin}}]}],

Sequence @@ sty, ptss}],
 FilterOptions[Graphics, opts], PlotRange->{{xmin,xmax}, {ymin,ymax}},
 Axes->Automatic, AxesOrigin->{xmin, ymin}];
If[showPts, {gg, pts1}, gg]]

PoincareSection[secondordereqn_, {t_, tmin_, tmax_}, 
     {x_, xmin_, xmax_}, {x_', ymin_, ymax_}, {x0_, xp0_}, opts___] :=
 ({sys, vbleChange} = ToSystem[
   secondordereqn, x[t], t, ShowAuxiliaryVariables -> True];
   PoincareSection[sys, {t, tmin, tmax}, {x, xmin, xmax},
   {x' /. vbleChange, ymin, ymax}, {x0, xp0}, opts])

ResidualPlot[x___] :=  (Message[VDS::insarg, ResidualPlot, 4]; 
    HoldForm[ResidualPlot[x]]) /; Length[{x}] < 4

ResidualPlot[x__] :=  (Message[VDS::nonopt, 4]; HoldForm[ResidualPlot[x]]) /; 

     Length[{x}] != 4 && Union[Head /@ Drop[{x}, 4]] =!= {Rule}

ResidualPlot[x__] :=  (Message[VDS::badarg]; 
            HoldForm[ResidualPlot[x]]) /; !ListQ[{x}[[4]]]

ResidualPlot[x__] :=  (Message[VDS::nonnum, {4,2}]; 
            HoldForm[ResidualPlot[x]]) /; !NumberQ[N[{x}[[4,2]]]]

ResidualPlot[x__] :=  (Message[VDS::nonnum, {4,3}]; 
            HoldForm[ResidualPlot[x]]) /; !NumberQ[N[{x}[[4,3]]]]
                       
Options[ResidualPlot] = {PlotType->Automatic};

ResidualPlot[eqn_Equal, soln_, unk_, {t_, tmin_, tmax_}, opts___] :=
  Module[{},
  halt = HoldForm[ResidualPlot[eqn,soln,unk,{t,tmin,tmax},opts]];
  If[FreeQ[eqn, Head[unk]'[t]] && FreeQ[eqn, Head[unk]''[t]], 
     Message[VDS::badform, t]; Return[halt]];
 ty = PlotType /. {opts} /. Options[ResidualPlot];
 If[!MemberQ[{Logarithmic, Automatic}, ty],
  Message[VDS::badopt, ty, PlotType]; Return[halt]];
 e = If[FreeQ[eqn, Derivative[2][Head[unk]][t]],
       Solve[eqn, D[unk, t]],
       Solve[eqn, D[unk, {t, 2}]]][[1,1]];
 Plot[Evaluate[If[ty === Logarithmic, Chop@Log[10, Abs[#]]&, Identity][
  e[[1]] - e[[2]]] /. Head[unk]->Function[t,soln]], {t, tmin, tmax}, 
  Evaluate[FilterOptions[Plot, opts]],
  PlotRange->All]]

           (* SecondOrderPlot *)

Options[SecondOrderPlot] = {PlotVariables->All};

SecondOrderPlot[e_?FreeOfSetQ, x___] :=  
  (Message[VDS::insarg, SecondOrderPlot, 3]; 
    HoldForm[SecondOrderPlot[e, x]]) /; Length[{x}] < 2

SecondOrderPlot[e_?FreeOfSetQ, x__] :=  (Message[VDS::badarg]; 
            HoldForm[SecondOrderPlot[e,x]]) /; !ListQ[{x}[[2]]]

SecondOrderPlot[e_?FreeOfSetQ, x__] :=  (Message[VDS::nonnum, {3, r[[1,
1]]+1}]; 
            HoldForm[SecondOrderPlot[e,x]]) /;
        (r=Position[(List[x][[2, {2,3}]]), 
            _?(!NumberQ[N@#]&), Heads->False])  != {{}}

(* Case for a PhasePlot *)
(* Case of no t-iterator *)
SecondOrderPlot[eqnn_?FreeOfSetQ /; Head[eqnn] =!= List, unk_, 
  {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, opts___Rule] := 
  SecondOrderPlot[eqnn, unk, {unk[[1]], 0, 1},  {x, xmin, xmax},
    {y, ymin, ymax}, opts]
 

SecondOrderPlot[eqnn_?FreeOfSetQ, unkks_, {t_, tmin_, tmax_}, 
 {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, opts___Rule] := 
  Module[{}, 
{eqn, unks} = If[!ListQ[eqnn], {{eqnn}, {unkks}}, {eqnn, unkks}];
halt = HoldForm[SecondOrderPlot[eqn, unks, {t, tmin, tmax}, 
    {x, xmin, xmax}, {y, ymin, ymax}, opts] ];
If[!ListQ[eqn], Message[VDS::badarg]; Return[halt]];
    
If[!CheckOptions[{SecondOrderPlot, PhasePlot, Graphics, NDSolve, 
  RungeKutta4, Euler, ContourPlot, Plot, DirectionFieldPlot,
  Equilibria, ColorParametricPlot, FlowParametricPlot,
  ProjectionPlot3D, Graphics3D, NullclinePlot}, opts], Return[halt]];
  
{sys, vbleChange} = ToSystem[eqn, unks, t, ShowAuxiliaryVariables -> True];
PhasePlot[Evaluate[sys],
    Evaluate[Join[unks, unks /. xx_[t] :> xx'[t] /. vbleChange ]],
    {t, tmin, tmax}, 
    Evaluate[Sequence @@ ({{x, xmin, xmax}, {y, ymin, ymax}} /. vbleChange)],
 Source -> SecondOrder,
 Evaluate[FilterOptions[PhasePlot, opts]],
 Evaluate[FilterOptions[ColorParametricPlot, opts]],
 Evaluate[FilterOptions[FlowParametricPlot, opts]],
 Evaluate[FilterOptions[ColorParametricPlot3D, opts]],
 Evaluate[FilterOptions[ProjectionPlot3D, opts]],
 Evaluate[FilterOptions[DirectionFieldPlot, opts]],
 Evaluate[FilterOptions[NDSolve, opts]],
 Evaluate[FilterOptions[ParametricPlot, opts]]]]
 

(* Main case: system or single equation *) 
SecondOrderPlot[eqnn_?FreeOfSetQ, unkks_, {t_, tmin_, tmax_}, opts___Rule] :=
 (
 
 halt = HoldForm[SecondOrderPlot[eqnn, unkks, {t, tmin, tmax}, opts]];
 
 If[!ListQ[eqnn] && Head[eqnn] =!= Equal, 
     Message[VDS::badarg]; Return[halt]];
 {eqn, unks} = If[!ListQ[eqnn] && Head[eqnn] === Equal, {{eqnn}, {unkks}},
   {eqnn, unkks}];
 If[!CheckOptions[{SecondOrderPlot, SystemSolutionPlot, Graphics, Plot, 
     NDSolve}, opts], Return[halt]];
plotvars = PlotVariables /. {opts} /. Options[SecondOrderPlot];
plotvars = plotvars /. (All | Automatic :> 
           {Join[Head /@ unks, Head[#]' & /@ unks]});

If[!ListQ[plotvars] || !(And @@ (ListQ /@ plotvars)), Message[
SecondOrderPlot::badvbl]; Return[halt]];
{sys, vbleChange} = ToSystem[eqn, unks, t, ShowAuxiliaryVariables -> True];
SystemSolutionPlot[sys,
  Evaluate[Join[unks, unks /. xx_[t] :> xx'[t] /. vbleChange ]],
    {t, tmin, tmax}, PlotVariables ->(plotvars /. vbleChange),
    Evaluate[FilterOptions[NDSolve, opts]],
    Evaluate[FilterOptions[Graphics, opts]],
    Evaluate[FilterOptions[SystemSolutionPlot, opts]]])


(* SystemSolutionPlot *)

SystemSolutionPlot[e_List?FreeOfSetQ, x___] :=  (Message[VDS::insarg, 
SystemSolutionPlot, 3]; 
    HoldForm[SystemSolutionPlot[e,x]]) /; Length[{x}] < 2

SystemSolutionPlot[e_List?FreeOfSetQ, x__] :=  (Message[VDS::nonopt, 3];
  HoldForm[SystemSolutionPlot[e,x]]) /; 
     Length[{x}] != 2 && Union[Head /@ Drop[{x}, 2]] =!= {Rule}

SystemSolutionPlot[e_List?FreeOfSetQ, x__] :=  (Message[VDS::badarg]; 
            HoldForm[SystemSolutionPlot[e,x]]) /; !ListQ[{x}[[1]]]

SystemSolutionPlot[e_List?FreeOfSetQ, x__] :=  (Message[VDS::badarg]; 
            HoldForm[SystemSolutionPlot[e,x]]) /; !ListQ[{x}[[2]]]

SystemSolutionPlot[e_List?FreeOfSetQ, x__] :=  (Message[VDS::nonnum, {3,2}]; 
           HoldForm[SystemSolutionPlot[e,x]]) /; !NumberQ[N[{x}[[2,2]]]]

SystemSolutionPlot[e_List?FreeOfSetQ, x__] :=  (Message[VDS::nonnum, {3,3}]; 
           HoldForm[SystemSolutionPlot[e,x]]) /; !NumberQ[N[{x}[[2,3]]]]

Options[SystemSolutionPlot] = {
 Background->GrayLevel[1],
 DefaultColor->Automatic,
 FastPlotting->False,
 InitialValues->{},
 PlotRanges->Automatic,
 PlotStyle -> {},
 PlotVariables->All,
 PlotLabels -> {},
 Rainbow->False,
 SaveLastPoint->False,
 SaveSolution->False,
 ShowNumberOfSteps -> False,
 SolutionName->"Solution",
 WorkingPrecision->$MachinePrecision};

SystemSolutionPlot[eqns_?FreeOfSetQ /; Length[eqns] > 1, unks_,
  {t_, tmin_, tmax_}, opts___] := Module[{},
halt = HoldForm@SystemSolutionPlot[eqns,unks,{t,tmin,tmax},opts];
  If[!CheckOptions[{SystemSolutionPlot, Graphics, Plot, NDSolve}, opts], 
   Return[halt]];

If[!ListQ[eqns], Message[VDS::badarg]; Return[halt]];
fixopts = Sequence @@ ({opts} /. RGBColor[a_, a_, a_] :> GrayLevel[a]);
    (* fixes DefaultColor bug in 2.2 *)

{showNumQ, pl, pRans, pSty, initVals, asp, saveSol, fastPlQ, defC,
 backCol, solName, plotvars, rain, wp, lastPtQ} = {ShowNumberOfSteps,
  PlotLabels, PlotRanges, PlotStyle, InitialValues, AspectRatio, SaveSolution,

  FastPlotting, DefaultColor, Background, SolutionName, PlotVariables,
  Rainbow, WorkingPrecision, SaveLastPoint} /. {fixopts} /.
  Options[SystemSolutionPlot];

If[!(MemberQ[{Automatic, All}, plotvars] || (
    ListQ[plotvars] && (And @@ (ListQ /@ plotvars)))),
   Message[VDS::badopt, plotvars, PlotVariables]; Return[halt]];
   
plotvars = plotvars /. All | Automatic -> {Head /@ unks};
If[{} != Complement[Flatten[plotvars], Head /@ unks],
   Message[SystemSolutionPlot::badvar]; Return[halt]];
nPlotVars = Length[plotvars[[1]]];
saveSol = saveSol || 
 ({opts} != {} && MemberQ[First /@ {opts}, SolutionName]);
defC = defC /. Automatic->(backCol /. z_?NumberQ :> 1. - z);

If[Length[plotvars] > 1 && !MemberQ[{Automatic, {}}, pSty],
    Message[SystemSolutionPlot::nopsty]; Return[halt]]; 
If[Length[plotvars] > 1 && rain,
    Message[SystemSolutionPlot::norain]; Return[halt]];
If[!StringQ[solName], Message[VisualDSolve::badname]; Return[halt]];
If[initVals=={}, Message[SystemSolutionPlot::noinit]; Return[halt]];
If[MemberQ[{All, Automatic}, pRans], pRans = pRans & /@ unks];
If[NumberQ[N[First[initVals]]], initVals = {initVals}];
initVals = N[initVals, wp];
num = Length[initVals];

If[!Apply[And, 
  Map[(Length[unks] <= Length[#] <= Length[unks] + 1)&, initVals]], 
  Message[SystemSolutionPlot::badln, Length[unks]]; Return[halt]];

eqns1 = NormalForm[eqns, unks];  (* to fix 2.2 bug *)

nSoln[iv_] := unks /. First[

If[Length[initVals] == 1, 
   NDSolve[Join[eqns1, (#[[1]] /. (t->First[iv])) == #[[2]] & /@ 
       Transpose[{unks, Rest@iv}]], unks, {t, tmin, tmax}, 
        FilterOptions[NDSolve, opts]], 
Check[stemp = NDSolve[Join[eqns1,
     (#[[1]] /. (t->First[iv])) == #[[2]] & /@ 
       Transpose[{unks, Rest@iv}]], unks, {t, tmin, tmax}, 
        FilterOptions[NDSolve, opts]],
 Message[VDS::mxstep, iv, t, getDomain[stemp[[1,1,2]]][[2]]]; 
 stemp, NDSolve::mxst]]];


solns = nSoln /@ (initVals /. 
     v_?(Length[#] == Length[unks] && !ListQ@First[#] &) :> Prepend[v, tmin]);


If[showNumQ, Print[
StringForm["The initial value `` required `` steps.", #[[2]],
   Length[GetPts[First[#[[1]]]]]]] & /@ Transpose[{solns, initVals}]];
pSty = fixSty[pSty, nPlotVars];

Clear[Evaluate[StringJoin["Global`", solName]]];
If[saveSol, Evaluate[ToExpression@StringJoin["Global`", solName]] =
   If[Length[initVals] == 1, First[solns], solns]];
If[lastPtQ, LastPoint = If[Length[initVals] == 1, First[solns], solns] /.
   t->N[tmax, wp]];
Clear[solNameTemp];
If[Length[#] == 1, First[#], #]& @ MapIndexed[(i = First[#2]; Show[
 If[fastPlQ,
    Table[Graphics[Map[Function[vb,
         (place = First@First@Position[unks, vb[t]];
     Append[If[i==1, 
          If[rain, 
            Prepend[pSty[[place]], Hue[1 - 0.8 j/num]], 
      pSty[[First@First@Position[#1,vb]]]], {}],
   Line[GetPtsFull[solns[[j, place]]]]])], #1], Axes->True], {j, num}],
 Table[Plot[Evaluate[solns[[j,First /@ (First /@ Position[unks, #[t]] & /@
 #1)]]],
{t, getDomain[solns[[j, 1]]][[1]], getDomain[solns[[j, 1]]][[2]]},
PlotLabel -> If[pl != {}, pl[[i]], None],
PlotStyle->If[i==1, If[rain, 
        Prepend[#, Hue[1 - 0.8 j/num]] & /@ pSty, pSty], {}], 
      DisplayFunction->Identity,
      Evaluate[FilterOptions[Plot, fixopts]]], {j, num}]],
  Evaluate[FilterOptions[Graphics, fixopts]],
  DisplayFunction->$DisplayFunction,
  PlotRange->pRans[[i]], DefaultColor->defC]) &, plotvars]];

Options[ToSystem] = {AuxiliaryVariables -> Automatic,
ShowAuxiliaryVariables -> False};

ToSystem[x___] :=  (Message[VDS::insarg, ToSystem, 3]; 
 HoldForm[VisualDSolve[x]]) /; Length[{x}] < 3
    
ToSystem[eqn_Equal, unk_, t_, opts___] := Module[{}, 
 {aux, auxQ} = {AuxiliaryVariables, ShowAuxiliaryVariables} /. {opts} /.
   Options[ToSystem] /. Automatic -> Unique[Head[unk]];
 rule = Head[unk]' -> aux; 
({{D[unk, t] == aux[t],
  aux'[t] == ((Solve[eqn, D[unk, {t, 2}]][[1,1,2]]) /.  -1. -> -1 /.
     rule)}, rule})[[If[auxQ, {1, 2}, 1]]]]

ToSystem[eqn_List, unk_List, t_, opts___] := (
 {aux, auxQ} = {AuxiliaryVariables, ShowAuxiliaryVariables} /. {opts} /.
   Options[ToSystem] /. Automatic -> (Unique[Head[#]] & /@ unk);
  rules = MapIndexed[Head[#1]' -> aux[[#2[[1]]]]&, unk];
 ({Join[MapIndexed[D[#1, t] == aux[[#2[[1]]]][t]&, unk],
        MapIndexed[aux[[#2[[1]]]]'[t] ==
    (Solve[#1, D[unk[[#2[[1]]]], {t, 2}]][[1,1,2]] /. rules)&, eqn]],
  rules})[[If[auxQ, {1, 2}, 1]]])

(* ****************** VisualDSolve  *************  *)

Options[VisualDSolve] =  {
  AspectRatio->1/GoldenRatio,
  Axes->True,
  Background->GrayLevel[1],
  ComputeWindow->Automatic,
  ContourLines->False,
  ContourShading->False,
  DefaultColor->Automatic,
  DirectionField->False,
  FastPlotting->False,
  FieldMeshSize->15,
  InflectionCurve->False,
  InflectionStyle->{AbsoluteDashing[{5, 5}]},
  InitialPointStyle->{GrayLevel[0], PointSize[0.021]},
  InitialValues->{},
  Isoclines->False,
  IsoclinePlotPoints->20,
  IsoclineShading->False,
  IsoclineStyle->{},
  IsoclineValues->{0},
  MaxBend->10,
  MaxSteps->500,
  Method->Automatic,
  PlotPoints->25,
  PlotStyle->{},
  PrecisionGoal->Automatic,
  Rainbow -> False,
  SaveLastPoint->False,
  SaveSolution->False,
  ShowInitialValues->False,
  ShowNumberOfSteps->False,
  SolutionName->"Solution",
  StayInWindow->False,
  SymbolicSolution->False,
  WindowShade->None,
  WorkingPrecision -> $MachinePrecision};

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((
    Message[VDS::insarg, VisualDSolve, 3]; 
    HoldForm[VisualDSolve[e,x]]) /; Length[{x}] < 2)

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::badarg]; 
            HoldForm[VisualDSolve[e,x]]) /; Head[e] =!= Equal)

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::badarg]; 
            HoldForm[VisualDSolve[e,x]]) /; !ListQ[{x}[[1]]])

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::badarg]; 
            HoldForm[VisualDSolve[e,x]]) /; !ListQ[{x}[[2]]])

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::nonnum, {2,2}]; 
           HoldForm[VisualDSolve[e,x]]) /; !NumberQ[N[{x}[[1,2]]]])

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::nonnum, {2,3}]; 
           HoldForm[VisualDSolve[e,x]]) /; !NumberQ[N[{x}[[1,3]]]])

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::nonnum, {3,2}]; 
           HoldForm[VisualDSolve[e,x]]) /; !NumberQ[N[{x}[[2,2]]]])

VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::nonnum, {3,3}]; 
           HoldForm[VisualDSolve[e,x]]) /; !NumberQ[N[{x}[[2,3]]]])
           
VisualDSolve[e_?FreeOfSetQ, x___]:=  ((Message[VDS::nonopt, 3]; 
      HoldForm[VisualDSolve[e,x ]]) /; 
     Length[{x}] != 2 && Union[Head /@ Drop[{x}, 2]] =!= {Rule})

VisualDSolve[equation_?FreeOfSetQ, {t_, tmin_, tmax_},
 {x_, xmin_, xmax_}, opts___] :=  
  Module[{nSoln},
halt =  HoldForm[VisualDSolve[equation,{t,tmin,tmax},{x,xmin,xmax},opts]];

If[!CheckOptions[{VisualDSolve, Graphics, RungeKutta4, Euler, NDSolve,
      ContourPlot, Plot, DirectionFieldPlot, FreehandAttempt}, opts], 
  Return[halt]];

(* the following line sets newOpts to be all the options that should
be used, whether they come from defaults or were added in via opts *)
newOpts = Sequence @@ Join[{opts}, Options[VisualDSolve]];

{showNumQ, fieldSteps, ptsQ, fieldQ, symQ, solSty, pts, isoQ, isoVals, infQ, 
ptCol,
 backCol, infSty, isoSty, asp, cppts, isoShade, saveSol, lastQ, fastPlQ, 
  maxS, method, defC, cLines, stayIn, solName, compWin, winShad, rain, wp} = 
   {ShowNumberOfSteps, FieldMeshSize, ShowInitialValues, DirectionField, 
SymbolicSolution,
    PlotStyle, InitialValues, Isoclines, IsoclineValues, InflectionCurve,
    InitialPointStyle, Background, InflectionStyle, IsoclineStyle,
   AspectRatio, IsoclinePlotPoints, IsoclineShading, SaveSolution,
   SaveLastPoint, FastPlotting, MaxSteps, Method, DefaultColor, ContourLines,
    StayInWindow, SolutionName, ComputeWindow, WindowShade, Rainbow,
    WorkingPrecision} /. 
   {opts} /. Options[VisualDSolve];
      
If[!fieldQ && pts == {} && !isoQ && !infQ, Message[VisualDSolve::nothng]; 
 Return[halt]];

Off[NDSolve::ndnum];
If[stayIn && symQ, Message[VisualDSolve::incons]; Return[halt]];
If[!StringQ[solName], Message[VisualDSolve::badname]; Return[halt]];
If[FreeQ[equation, x'[t]], Message[VDS::badform, t]; Return[halt]];

eqn = NormalForm[equation, x[t]];
saveSol = saveSol || MemberQ[First /@ {opts}, SolutionName];
If[saveSol && symQ, Message[VisualDSolve::incons1]; Return[halt]];
If[symQ && lastQ, Message[VisualDSolve::incons2]; Return[halt]];
  
tminn = Min[tmin, tmax]; tmaxx = Max[tmin, tmax];
xmaxx = Max[xmin, xmax]; xminn = Min[xmin, xmax];

{cwx0, cwx1} = N[compWin /. Automatic -> {xminn, xmaxx}];
{tR, xR} = {tmaxx - tminn, xmaxx - xminn};
axLab = AxesLabel /. {opts} /. AxesLabel->None;
pRange = PlotRange /. {opts} /. PlotRange->{{tminn, tmaxx}, {xminn, xmaxx}};
If[Depth[pts] == 2 && Head[pts] =!= Grid, pts = {pts}];

pts = pts /. Grid[solnMesh_] :> ({tstep, xstep} = N[{tR, xR}/(solnMesh + 1)];
      Flatten[Table[
         {tt, xx, {tminn, tmaxx}}, {tt, tminn+tstep, tmaxx-tstep, tstep}, 
                   {xx, xminn+xstep, xmaxx-xstep, xstep}], 1]);

asp = asp /. Automatic->xR/tR;
transformed = False;
method = method /. Automatic->NDSolve;
solSty = fixSty[solSty, Length[pts]];
If[!ListQ[ptCol], ptCol = {ptCol}];
If[method =!= NDSolve && stayIn, Message[VisualDSolve::bdmet]; Return[halt]];

stopOpt = If[stayIn && ver3Q, 
    StoppingTest -> !(cwx0 <= x[t] <= cwx1), WorkingPrecision -> wp];

nSoln[{tt_, xx_, {t1_, t2_}}] := (
 If[stayIn && {t1, t2} == {N[tmin, wp], N[tmax, wp]} && !transformed &&
  !ver3Q, 
  transformed = True;
  eqn = (eqn[[1]] == If[cwx0 <= x[t] <= cwx1, Evaluate[eqn[[2]]]]) ];

If[Length[pts] == 1, 
      x[t] /. method[{eqn, x[tt] == xx}, x[t], {t, t1, t2}, 
      Evaluate@FilterOptions[method, newOpts], Evaluate[stopOpt]],

      x[t] /. Check[stemp = method[{eqn, x[tt] == xx}, x[t], {t, t1, t2}, 
      Evaluate@FilterOptions[method, newOpts], Evaluate[stopOpt]],
      Message[VDS::mxstep, {tt, xx}, t, getDomain[stemp[[1,1,2]]][[2]]];
        stemp,
      NDSolve::mxst]]);

defC = defC /. (Automatic->backCol /. z_?NumberQ :> 1. - z);
  (* Because of a bug in DefaultColor, this is programmed
  manually to be complementary to the background. *)

{fieldSteps, solnMesh} = 
   If[NumberQ[#], {#, #}, #]& /@ {fieldSteps, solnMesh};
   
If[TensorRank[pts] == 1, pts = {pts}];
pts = N[pts /. {u__?(NumberQ[N[#]]&), 
                  {t1_?(NumberQ[N[#]]&), t2_?(NumberQ[N[#]]&)}} :> 
     {u, {t1, t2, xx}} /.
     {u__?(NumberQ[N[#]]&)} :> {u, {tmin, tmax}} /.
     {u__, {t1_, t2_, xx}} :> {u, {t1, t2}}];
If[TensorRank[pts] == 3, pts = pts[[1]]];

(* Now we attempt to find a symbolic solution, if so requested *)
If[pts == {{}}, symQ = False];
symOut = {};
symTry =  If[symQ,  DSolve[eqn, x[t], t]];
symWorked = If[Length[symTry] > 1, 
  Message[VisualDSolve::manycs]; symOut = symTry; False,
   (symQ && FreeQ[symTry, Integrate] && FreeQ[symTry, DSolve] &&
     (constantSoln = Check[First /@  Solve[ (
     x[t] /. First @ First @ symTry /. t->#[[1]]) == #[[2]] , C[1]] & /@
     pts, "Solve failed"];
   (!StringQ[constantSoln]) && Union[Length /@ constantSoln] == {1}))];

If[symWorked, symOut = symTry];
If[symQ && symOut == {}, Message[VisualDSolve::nosym]];

If[symWorked && Length[symTry] != 1,
  Message[VisualDSolve::manycs, symTry]; symWorked = False];

Clear[Evaluate@StringJoin["Global`", solName]];
If[saveSol, solNameTemp = {}];
If[stayIn, WindowExitData = {}];
If[lastQ, LastPoint = {}];

If[symOut != {}, {#, symOut /. C[1]->C}, 
  If[stayIn, {#, WindowExitData}, #]] & @ Show[
If[winShad === None, {}, 
Graphics[{{winShad, Rectangle[{tminn, xminn}, {tmaxx, xmaxx}]}}]],
   
If[isoShade, Graphics[ContourPlot[Last[eqn] /. x[t]->x,
       {t, tminn, tmaxx}, {x, xminn, xmaxx},
	    PlotPoints->cppts,
	    ContourShading->True, Frame->False, PlotRange->All, 
	    ColorFunction->(GrayLevel[# .6 + .4] &),
     	DisplayFunction->Identity,
     	Evaluate[FilterOptions[ContourPlot, newOpts]],
     Contours->isoVals ]] /. 
     	   If[(ContourLines /. {newOpts}) === True, {}, _Line->{}] , {}], 
 If[fieldQ, N[DirectionFieldPlot[Last[eqn] /. x[t]->x, 
  {t, tminn, tmaxx, tR/fieldSteps[[1]]}, 
  {x, xminn, xmaxx, xR/fieldSteps[[2]]}, 
   Evaluate[Sequence[FilterOptions[DirectionFieldPlot, newOpts], 
                     FilterOptions[Graphics, newOpts]]],
   DisplayFunction->Identity, Source->SingleEqn, 
   PlotRange->pRange]], {}],

	 If[isoQ, Graphics[
	   ContourPlot[Last[eqn] /. x[t]->x,
       {t, tminn, tmaxx}, {x, xminn, xmaxx},
	    PlotPoints->cppts, PlotRange->All,
	    ContourLines->True,
	    Contours->isoVals, ContourStyle->isoSty,
	    ContourShading->False,
	    DisplayFunction->Identity,
	    Evaluate[FilterOptions[ContourPlot, newOpts]]]], {}],   
	    
	    Graphics[{If[infQ, Graphics[ContourPlot[ 
 Evaluate[D[Last[eqn], t] /. ToRules[eqn] /. x[t]->x],
    {t, tminn, tmaxx}, {x, xminn, xmaxx},
      Evaluate[FilterOptions[ContourPlot, newOpts]],
      ContourStyle->infSty,
       Contours->{0}, DisplayFunction->Identity]][[1]], {}]}],

If[pts == {{}}, {},
 If[symWorked, MapIndexed[
    Plot[#1[[1]], {t, #1[[2, 1]], #1[[2, 2]]},
    PlotPoints->cppts,
    PlotStyle->Evaluate[If[rain,
       {Prepend[solSty[[First[#2]]], Hue[1 - 0.8 #2[[1]] / Length[pts] ]]},
        solSty[[#2]] ] ],
    Evaluate[FilterOptions[Plot, newOpts]], 
    DisplayFunction->Identity]& , 
       Transpose[{x[t] /. symTry /. constantSoln, Last /@ pts}]], 
  MapIndexed[(sol = First[nSoln[#1]]; 
  If[method === NDSolve && showNumQ, Print[StringForm[
   "The initial value `` required `` steps.", 
      Drop[#1, -1], Length[GetPts[sol]]]]];
  
  If[saveSol, AppendTo[solNameTemp, sol]];
  If[stayIn, 
    AppendTo[WindowExitData, {#1[[1]], #1[[2]], getDomain[sol]}]];
  If[lastQ, AppendTo[LastPoint, sol /. t->Min[
  getDomain[sol][[2]], #1[[-1, 2]]]]];
  If[! ((sol /. InterpolatingFunction[___][t] :> True) &&
   getDomain[sol][[2]] > getDomain[sol][[1]]), 
       Message[VisualDSolve::failure]; Graphics[{}],
  If[fastPlQ, Graphics[{Append[
    If[rain, Prepend[solSty[[First[#2]]], Hue[1 - 0.8 #2[[1]] / Length[pts]
 ]],
      solSty[[First[#2]]] ],
  Line[GetPtsFull[sol]  ]]}],
  Plot[sol, 
   {t, Max[getDomain[sol][[1]], #1[[-1, 1]]], 
       Min[getDomain[sol][[2]], #1[[-1, 2]]]},
   PlotStyle->Evaluate[If[rain,
       {Prepend[solSty[[First[#2]]], Hue[1 - 0.8 #2[[1]] / Length[pts] ]]},
        solSty[[#2]] ] ],
   Evaluate[FilterOptions[Plot, newOpts]], 
   DisplayFunction->Identity ]]]) &,   pts]   ]], 
On[NDSolve::ndnum];
If[Length[solNameTemp] == 1, solNameTemp = First[solNameTemp]];
Evaluate[ToExpression@StringJoin["Global`", solName]] = solNameTemp;
Clear[solNameTemp];

If[Length[LastPoint] == 1, LastPoint = First[LastPoint]];
 Graphics[{{defC, AbsoluteThickness[1], 
    Line[{{tminn, xmaxx}, {tmaxx, xmaxx}, {tmaxx, xminn}, {tminn, xminn}, 
              {tminn, xmaxx}}]},
    If[ptsQ, Append[ptCol, Point /@ (Drop[#, -1] & /@ pts)], {}]}],
 AxesStyle->backCol,
 Frame->False,
 DefaultColor->defC,
 PlotRange->pRange,
 AxesLabel->axLab,
 Evaluate@FilterOptions[Graphics, newOpts],
 Ticks->Map[Flatten, 
     FullOptions[Graphics[{Line[{{tminn, xminn}, {tmaxx, xmaxx}}]},
     Axes->Automatic, AxesOrigin->{tminn, xminn},
					AspectRatio->asp], Ticks]  /. {a_, b_, _, _} :> If[b == "", {}, a] ],
 AxesOrigin->{tminn, xminn},
 DisplayFunction->$DisplayFunction]]

End[]; EndPackage[];