% This is a LaTeX file.  To use it you must first process it with
% LaTeX. Then you can render it with whatever printer you want.  If
% you don't have TeX, then you should.  It is free and anyone doing
% mathematics should have it.
\documentstyle{article}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\beqa}{\begin{eqnarray}}
\newcommand{\eeqa}{\end{eqnarray}}
\newcommand{\bvb}{\begin{verbatim}}
\newcommand{\evb}{\end{verbatim}}
\newcommand{\beqann}{\begin{eqnarray*}}
\newcommand{\eeqann}{\end{eqnarray*}}
\newcommand{\nn}{\mbox{${\nonumber}$}}
\newcommand{\labeq}[1]{\label{eq:#1}}
\newcommand{\refeq}[1]{(\ref{eq:#1})}
\newcommand{\tc}[1]{\addcontentsline{toc}{subsection}{#1}}
\newcommand{\TCC}[1]{\medskip

\addcontentsline{toc}{subsubsection}{#1}
\noindent{\bf {#1}.}}
	
\newcommand{\tcc}[1]{\addcontentsline{toc}{subsubsection}{#1}}
\title{XPPAUT3.05 -- the differential equations tool}
\author{Bard Ermentrout}
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Introduction} XPP (XPPAUT is another name; I will use the two
interchangeably)
  is a tool for solving differential
equations,
 difference 
equations, delay equations, functional equations, boundary value
problems, and stochastic equations. 
 It evolved 
from a chapter written by John Rinzel and myself on the qualitative
 theory of nerve membranes and eventually became a commercial product
 for MSDOS computers called PHASEPLANE.  It is now available as a
program running
 under
 X11 and UNIX.  

 The code 
brings together a number of useful algorithms and is extremely portable. 
 All the graphics and interface are written completely in Xlib which explains 
the somewhat idiosyncratic and primitive widgets interface.  

XPP contains the code for the popular bifurcation
program, AUTO.  Thus, you can switch back and forth between XPP and
AUTO, using the values of one program in the other and vice-versa.  I
have put a ``friendly'' face on AUTO as well.  You do not need to know
much about it to play around with it.

XPP has the capabilities for handling up to 300 differential equations. 
 There are solvers for delay and stiff differential equations a well as
 some code for boundary value problems.  Difference equations are also
 handled.  Up to 10 graphics windows can be visible at once and a variety
 of color combinations is supported.  PostScript output
 is supported.  Post processing is easy and includes the ability to make
 histograms, FFTs and applying functions to columns of your data.
Equilibria
 and linear stability as well as one-dimensional invariant sets can be
 computed.  Nullclines and flow fields aid in the qualitative understanding
 of two-dimensional models.  Poincare maps and equations on cylinders and
 tori are also supported.  Some useful averaging theory tricks and various
 methods for dealing with coupled oscillators are included primarily because
 that is what I do for a living. Equations with Dirac delta functions
are allowable.

With Version 3.0, I have added an animation package that allows you to
create animated versions of your simulations, such as a little
pendulum moving back and forth or lamprey swimming. The new window for
animations is produced by invoking the {\tt (V)iew axes} {\tt (T)oon}
menu item.  See section \ref{toon} for complete info.

 
I will assume that you are well versed in the theory of ordinary differential
 equations although you need not be to use the program.  
There are a number of useful features designed for people who use
dynamical systems to {\em model} their experiments.  There is a
curve-fitter based on the Marquardt-Levenberg algorithm which lets you
fit data points to the solutions to dynamical systems.  
Gnuplot-like graphics and
support for some graphics objects such as text, arrows, and
pointers are part of the package. You can also import bifurcation
curves as part of your graphs.
 It is possible to automatically generate ``movies'' of
three-dimensional views of attractors or parametric changes in the
attractor as some parameters vary.    I have
also included a small preprocessing utility that allows one to create
files for large systems of coupled equations. 

There are a number of other such programs available, but they all seem to
 require that your problems be compiled before using them.  XPP does not;
 I have devised a simple and fairly fast formula compiler that is based
 on the idea of the inner interpretor used in the language FORTH (which
 remains my first love as far as language is concerned)  Fear not, the
 differential equations and boundary conditions and other formulae are
 written in usual algebraic notation.  However, in order to run big
problems very quickly, I have written the code so that it is possible
to create a library that can be linked to your problem and thus create
a binary with the right-hand sides compiled.  This can run much faster
than the parsed code.  See the notes on this at the end of the
document. 




XPP has been successfully compiled on a SPARC II under OpenLook, a SPARC
 1.5 running generic X, a NeXT running X11R4, a DEC 5000, a PC using
Linux, and SGI and an HP 730.
  I cannot vouch for other platforms but it has been compiled on the
 IBM RS6000. Building XPP requires only the standard C compiler,
 and Xlib.  Look at the any README files that come with the
distribution for solutions to common compilation problems.



The basic unit for XPP is a single ASCII file (hereafter called an ODE
file) that has the equations, parameters,
 variables, boundary conditions, and functions for your model. You can
also include numerical parameters such as time step size and method of
integration although these can also be changed within the program.
The graphics and postprocessing are all done within
 the program using the mouse and various menus and buttons.  The impatient
 user should look at some sample {\tt *.ode} files instead of actually reading 
the documentation.  In addition, XPP uses a file called {\tt default.opt} that
 describes initialization and memory options for the program. This is
not necessary as all of the information contained in this file can now
be included in ODE file.
\addcontentsline{toc}{subsection}{Interface}
\begin{center}
{\Large Notes on the Interface }
\end{center}
This interface is crude by most standards and it is rather ugly as
well. I had originally thought to do the whole thing using the
OpenLook Widget Set but found that many of the systems to which I
wanted to port XPP did not have the set.  The present design is
completely portable.  Someday I will redo it in something like Motif.   

The text editing is also somewhat restricted.  There is as yet no cut
and paste.  However,
one can use the {\tt Home,End} keys to move the text cursor to the
beginning and end of the line.  The left and right arrow keys let you
move back and forth.  You can  insert or delete text at any point.
{\bf NOTE:} Be careful with the BackSpace and Delete keys as onj some
systems they are mapped differently.  One of them causes the whole
line to be erased and the other just erases a single character.
Almost every command has a keyboard shortcut.  These are given below.


\tc{Disclaimer}
\begin{center} 

\Large Disclaimer \\

\end{center}

XPP is distributed as is.  The author makes no claims as to the performance
 of the program.  Anyone is allowed to modify and distribute XPP as long as
 the original code is also made available.  

\tc{Acknowledgments}
\begin{center}
Acknowledgements
\end{center}
The porting of XPP to UNIX and X has been partially supported by NSF
in the sense that they gave me the SPARC 2.  MSRI in Berkeley gave me
use of a lovely office facing the bay and a workstation on which the
majority of the port was made in the Summer of 1992.  The original
version of PhasePlane benefitted by countless colleagues whose endless
requests for features and whose constant use led to the version that
you see now.  Of these I would like to specifically cite John Rinzel and
Artie Sherman at the NIH for being the main guinea pigs for the
earliest versions.

Please let me know of any bugs or other stuff that you'd like to see
incorporated into XPP.  I will usually fix them quickly.  
  
My EMAIL address is bard@mthbard.math.pitt.edu.
\tc{World Wide Web Tutorial}
\begin{center}
{\bf NOTE}
\end{center}
The easiest way to get a thorough understanding of the program as well
as a short tutorial in dynamical systems is to use the World Wide Web
tutorial which can be accessed from my home page at 
http://www.pitt.edu/$\equiv$phase. This tutorial is geared toward
computational neuroscientists (in the choice of problems) but provides
a fairly detailed introduction to the program.  


\section{ODE Files}

{\bf NOTE.}  {\em The latest version of XPP has a greatly improved
parser. Any versions 1.6 and above will have this parser built into
them.  The new parser is compatible with the old format for ODE files
so that it will read any old style files.  Thus, for users of previous
versions, you can still use all of the old files.
However, the old style parser will no longer be described in the
documentation.}  A command line option lets you convert old-style to
new style format.


ODE files are ASCII readable files that the XPP parser reads to create
 machine usable code. Lines can be continued with the standard
backslash character, however, the total length of any line cannot
exceed 256 characters.

\bigskip
\tc{Example}
{\bf Example.} 
I will start with a very simple example to get you up and running. The
model is the periodically driven Fitzhugh-Nagumo equation:
\beqa
dv/dt &=& f(v)-w+s(t)+I_0 \\
dw/dt &=& \epsilon(v-\gamma w) \\
f(v) &=& v(1-v)(v-a) \\
s(t) &=& \alpha \sin \omega t
\eeqa

Here is the ODE file:
\begin{verbatim}
# Forced Fitzhugh-Nagumo fhn.ode 
dv/dt = f(v)-w+s(t)+I_0
dw/dt = eps*(v-gamma*w)
f(v)=v*(1-v)*(v-a)
s(t)=al*sin(omega*t)
param a=.25,eps=.05,gamma=1,I_0=.25
param al=0,omega=2
@ total=100,dt=.2,xhi=100
done
\end{verbatim}

The file is pretty self-explanatory.  The first line cannot contain a
number as its first character (this makes the parser think that the
format of the ODE file is the old style.)  The last line should be the
word ``done.''  The names of all parameters must be declared with
optional values (the default sets them to 0.) There can be as many as
you can fit on each line (up to 200 parameters) but they must be
separated by commas or spaces and the ``='' sign must have {\em no} spaces on
either side of it.  

You could optionally include initial data by adding either of the
following sets of lines to the file:
\begin{verbatim}
init v=.25,w=.3
\end{verbatim}
or
\begin{verbatim}
v(0)=.25
w(0)=.3
\end{verbatim}

As you have probably guessed, comments have the form:
\begin{verbatim}
# This is a comment
\end{verbatim}

The ``@'' sign tells XPP that you want to preset some of the internal
parameters for numerical integration and graphing.
  In this case, we have told XPP
to integrate the equations until t=100 with a timestep of 0.25 and to
set the high value for the x-axis to 100.  You can, of course, change
all these internal options from within XPP; this provides an easy way
``set'' up the problem for ``one-button'' operation. 

\tc{Quick exploration.}
\begin{center} {\Large Quick exploration} \end{center}

Once you have written an ODE file, you can run XPP by typing
\begin{verbatim}
xpp ffhn.ode
\end{verbatim}
where {\tt ffhn.ode} is the filename you created.  You may have to
type
\begin{verbatim}
xppaut ffhn.ode
\end{verbatim}
depending on what you named the program. (The distributed source
produces a binary called {\tt xppaut} , but I usually create a shell
script called {\tt xpp} that contains the following line
\begin{verbatim}
xppaut $1 -xorfix
\end{verbatim}
thus alleviating problems my display has with the exclusive or drawing
required for ``zooming'' in and other things.)

When you type one of the above lines, then the program fires up. It
reads the ODE file and if there are any errors, it reports them and
returns to the command line.  Assuming you haven't written the ODE
file wrong (which you shouldn't have since it is included in the
source) you should get the main window.  Depending on your X windows
system, the other windows may or may not be iconified.  Almost all
command have keyboard shortcuts which are either the first letter of
the command or the letter in parentheses or capitalized. Thus, you can
either use the mouse to click on the command or you can use the
keyboard to choose the command.  To solve the differential equation
with the current parameters, click on {\tt Initialconds} and then {\tt
Go}. You will see the variable $V(t)$ plotted across the screen as a
function of time.  (Instead of using the mouse, you could type {\tt I
G} as a keyboard shortcut.)  Click on {\tt Xivst}.  When the prompt
comes up backspace over {\tt V}, and type in {\tt w} and {\tt Enter.}
The variable $w$ will be plotted versus time.  Note that the vertical
axis of the window is automatically adjusted.  Click on {\tt Viewaxes}
to choose the view and when the new window pops up fill it in as follows:
\begin{itemize} {\tt
\item X-axis: V
\item Y-axis: W
\item Xmin: -.5
\item Ymin: 0
\item Xmax: 1.5
\item Ymax: 1
\item Xlabel: V
\item Ylabel: w
}
\end{itemize}
and then click on {\tt Ok.}  The phase-plane will be drawn showing a
limit cycle. Click {\tt Nullclines} and then {\tt New} to draw the
nullclines. Click {\tt Text,etc} then {\tt Text} and type in
``V-nullcline'' followed by {\tt Enter} at the prompt.  Accept the
defaults for text size and font by typing {\tt Enter} twice.  Move the
mouse pointer to the cubic-like curve and click the button.  The text
should appear on the screen.  Repeat this but type in ``w-nullcline''
for the text.  Click {\tt Graphic stuff} and then {\tt Postscript}.
Accept the defaults and a hardcopy postscript file will be produced
which you can view or printout on an appropriate printer.  Click {\tt
File} and then {\tt Quit} and answer {\tt Yes} to exit XPP.  

A much more extensive tutorial is available on the World Wide Web (see
above). This document is mainly a reference to all the features (bugs
:) of XPP.

 
\begin{center}
{\Large ODE File format }
\end{center}
\tc{ODE File format} 
ODE files consist of ascii readable text which XPP uses to describe
the program it wants to solve.  The line length is limited to 256
characters total. Individual lines can be continued with the UNIX
backslash character, $\backslash.$ ODE files have any combination of the
following lines.  The order is not too important but can matter (see
below).
 

\begin{verbatim}
# comment line - name of file, etc   
options <filename>
...
d<name>/dt=<formula>
<name>'=<formula>
...
<name>(t)=<formula>
...
volt <name>=<formula>
...
<name>(t+1)=<formula>
...
markov <name> <nstates>
{t01} {t02} ... {t0k-1}
{t10} ...
...
{tk-1,0} ...
...
aux <name>=<formula>
...
<name>=<formula>
...
parameter <name1>=<value1>,<name2>=<value2>, ...
...
wiener <name1>, <name2>, ...
...
number <name1>=<value1>,<name2>=<value2>, ...
...
<name>(<x1>,<x2>,...,<xn>)=<formula>
...
table <name> <filename>
...
table <name> % <npts> <xlo> <xhi> <function(t)>
...
global sign {condition} {name1=form1;...}
...
init <name>=<value>,...
...
<name>(0)=<value>
...
bdry <expression>
...
command[i1..i2] ...
...
name[i1..i2] ...
...
# comments
...
@ <name>=<value>, ...
...
set <name> {x1=z1,x2=z2,...,}
..
done
\end{verbatim}
       



The typical ODE file contains some or all of the above types of lines.
Continuous variables, auxiliary quantities, and Markov variables are
all plottable quantities in XPP.  That is, once you have solved your
equation, you can plot or view any of the continuous and Markov
variables or the auxiliary quantities. 

\TCC{Options files}
XPP uses a bunch of defaults when it is started up by looking for a
file called ``default.opt.'' If it cannot find it, it uses internal
options.  Alternatively, you can tell XPP the name of the options file
you want to use.  The description of these files is below.
The format for such a statement is:
\bvb
option <filename>
\end{verbatim}
which loads the options file specified in {\tt<filename>}. You will
probably not want to use this very much as you can now specify all of
the parameters in the options file within your ODE file by using the
``@'' symbol. 

\TCC{Defining Continuous Variables}
Variables are the quantities you wish to integrate in time. There are
two types of variables: (i) continuous and (ii) Markov. I will first
describe continuous variables.  Variable names (as
can all names in XPP) can have up to 9 letters each.  XPP is case
insensitive. Any combination of letters and numbers is valid as is the
underscore, ``\_''. There are 5 ways that you can tell XPP the names of
the continuous variables and their right-hand sides. The following three are
equivalent:
\begin{verbatim}
d<name>/dt=<formula>
<name>'=<formula>
<name>(t)=<formula>
\end{verbatim}
Here {\tt <name>} is the name of the variable.  The last version is
for notational convenience only sincd the {\tt dx/dt} notation makes
no sense for discrete dynamical systems.  The {\tt <formula> } is
exactly that, the formula for the right-hand sides of the equations.
These equations can appear anywhere in your file and in any order.
However, the order in which they are written determines the order in
which they appear in the Data Browser (see below.) 

The fourth way of defining a continuous variable is:
\begin{verbatim}
<name>(t) = <formula>
\end{verbatim}
which tells XPP that this defines a Volterra integral equation. (It is
distinguished from a definition of some function of the dummy variable
{\tt t} by the presence of an integral operator ({\tt int\{} or {\tt
int[} ) in the right-hand side. 
(see below).
For example, the convolution equation:
\[
 v(t) = \exp(-t) + \int_0^t e^{-(t-s)^2}v(s) ds
\]
would be written as:
\begin{verbatim}
v(t) = exp(-t) + int{exp(-t^2)#v}
\end{verbatim}
Integro-differential equations use the {\tt dv/dt} etc notation so
that
\[
 \frac{dv(t)}{dt} = -v(t)+ \int_0^t e^{-(t-s)^2}v(s) ds
\]
becomes
\begin{verbatim}
dv/dt= -v+int{exp(-t^2)#v}
\end{verbatim}
In the event that the right-hand side does not contain any integral
operator (as would be the case, for example, if there was a fixed or
hidden variable definition) then you can force the parser into making
the equation a Volterra integral equation by typing
\begin{verbatim}
volt v= exp(-t)+int{exp(-t^2)#v}
\end{verbatim}
{\bf NOTE.} In this format you do not write {\tt v(t)=...} but just
{\tt v=...}

\TCC{Pseudo-arrays}
XPP gives you the option of defining many variables at once using an
array-like declaration which XPP expands into a set of
declarations. For example, you could declare 10 equations with the
following command:
\begin{verbatim}
x[1..10]'=-x[j]
\end{verbatim}
and XPP would internally expand this is
\begin{verbatim}
dX1/dt=-X1
dX2/dt=-X2
dX3/dt=-X3
dX4/dt=-X4
dX5/dt=-X5
dX6/dt=-X6
dX7/dt=-X7
dX8/dt=-X8
dX9/dt=-X9
dX10/dt=-X10
\end{verbatim}
Thus, you can make networks and discretizations of PDE's compactly.
Note the appearance of the expression {\tt [j] }.   This is expanded
by XPP to the value of the index. Similarly, the following are
allowable indices:
\begin{verbatim}
[j+n]
[j-n]
[j*n]
\end{verbatim}
where {\tt n } is any integer.  The {\tt [1..10] } notation tells XPP
to start {\tt j} at 1 and go to 10.  You can start with any
nonnegative integer and end with any. The first number can be less
than or greater than the second. XPP does not treat arrays in
any efficient manner, it is as if you defined 10 or whatever
variables. The names of the variables are the root name with the index
appended.

\TCC{Markov Variables}  
Markov variables are finite state quantities that randomly flip from
one integer state to another according to the transition probability table
that is given in the ODE file. They are treated like variables in that
they have initial conditions and are accessible to the user. 
Markov variables are declared as
\begin{verbatim}
markov <name> <nstates>
{t01} {t02} ... {t0k-1}
{t10} ...
...
{tk-1,0} ...
...
\end{verbatim}
Each Markov variable has
its own line which must begin with the letter ``m''. (All other
letters are ignored, but for readability, it is best to write it out.)
  The name of the variable, {\tt <name>}, and   the
number of states, {\tt <nstates>}, are included on the first line.  
The possible values are 0,1, and so
on up to $k-1$ where $k$ is the number of states.

A
finite state Markov variable with $k$ states must have associated with
it a $k\times k$ matrix, $T_{ij}$ which contains the probability of going from
state $i$ to state $j$ per unit of time.  Thus the effective
transition probability is {\tt DeltaT} times $T_{ij}$; the
larger is {\tt DeltaT} the higher the probability.  Since the
probabilities must add to 1 in any row, the diagonal terms are
automatic and ignored by XPP. 
 If there are $k$ states to the variable, then there must be
$k$ rows following the declaration of the transition matrix.  Each row
contains $k$ entries delimited by curly brackets and separated by
spaces. For example, suppose $z$ is a two-state variable with
transition probabilities, $P_{ij}$ then it would be defined by:
\begin{verbatim}
markov z 2
{0} {P01}
{P10} {0}
\end{verbatim}
where {\tt P01, P10} are any {\em algebraic expression or number involving
the parameters and variables} of the XPP file. Note that this implies
that the state transitions can be dependent on any other quantities.

At each output time step, the probabilities are computed, multiplied by
the timestep, and a random number is chosen.  If it is in the appropriate
range, then the transition will be made.  Transitions of several such
variables are made in parallel and then each is updated. 

\TCC{Auxiliary quantities} In many cases, you might one to keep track
of some combination of your variables.  For example, you might want to
track the potential energy of a damped pendulum as it swings. They are
declared as:
\begin{verbatim}
aux <name>=formula
\end{verbatim}
where {\tt <name>} is the name of the quantity and {\tt <formula>} is
the formula for it.  {\em Note that a formula cannot refer to an
auxiliary named quantity; use fixed or hidden variables for this.}
An example using the auxiliary quantity is the damped pendulum:
\[
ml\frac{d^2x}{dt^2} = -mg\sin x -\mu\frac{dx}{dt}
\]
with potential energy:
\[
P.E. = mg(1-\cos x)
\]    
and kinetic energy
\[
K.E.= \frac{1}{2}ml(\frac{dx}{dt})^2.
\]
Since XPP solves systems of first order equations, this is first
converted and results in the ODE file:
\begin{verbatim}
# damped pendulum 
dx/dt = xp
dxp/dt = (-mu*xp-m*g*sin(x))/(m*l)
aux P.E.=m*g*(1-cos(x))
aux K.E.=.5*m*l*xp^2
param m=10,mu=.1,g=9.8,l=1
done
\end{verbatim}
where I have also given some values to the parameters. 

As with differential equations, you can also define many auxialiar
variables at once with a statement like:
\begin{verbatim}
aux r[1..10]=sqrt(x[j]^2+y[j]^2)
\end{verbatim}
which will be expanded in the obvious fashion.
 
\TCC{Hidden or fixed quantities} XPP allows you to define intermediate
quantities that can be used in the right-hand sides of the
equations. They are kept internally by XPP and here, the order in
which they are declared matters.  They are evaluated in the order in
which they are defined, so earlier defined ones should not refer to
later defined ones.  The format is:
\begin{verbatim}
<name> = <formula>
\end{verbatim}
They are most useful if you want to use a complicated quantity in
several right-hand sides.  The {\tt <name>} is kept internal to XPP
and their values are not stored (unlike variables). {\em Note that
they are different from functions which can take arguments and are not
hidden from the user.} 

For example, in the pendulum model above, you might also want the
total energy:
\[
T.E. = K.E. + P.E.
\]
Now, as I remarked above, yopu cannot just add another auxiliary
variable using {\tt P.E.} and {\tt K.E.} since they are not known to
XPP. But why compute them twice.  Here is how to use fixed variables
in this example.

\begin{verbatim}
# damped pendulum pend.ode
dx/dt = xp
dxp/dt = (-mu*xp-m*g*sin(x))/(m*l)
pe=m*g*(1-cos(x))
ke=.5*m*l*xp^2
aux P.E.=pe
aux K.E.=ke
aux T.E=pe+ke
param m=10,mu=.1,g=9.8,l=1
done
\end{verbatim}
Both energies are only computed once. (For this example, the
performance difference for computing the additional quantity is
negligible, but for more complex formulae, fixed quantities are useful.) 

Hidden variables can also be declared in groups like ODEs:
\begin{verbatim}
ica[1..10]=gca*minf(v[j])*(v[j]-eca)
\end{verbatim}
this is expanded into 10 declarations:
\begin{verbatim}
ica1=gca*minf(v1)*(v1-eca)
ica2=gca*minf(v2)*(v2-eca)
...
\end{verbatim}


\TCC{Parameters, Wiener parameters, Numbers} 
Parameters are named quantities that represent constants in your ODE
and which you can change from within the program.  The format is:
\begin{verbatim}
parameter <name1>=<value1>, <name2>=<value2>,...
\end{verbatim}
There can be many declarations on each line.  It is very important
that there be {\em no spaces} between the {\tt <name>} the {\tt =}
sign, and the {\tt value} of the parameter.  Without an {\tt =} sign
and a value, the parameter is set to zero by default.  

Numbers are just like parameters but they are ``hidden'' from the
user; they do not appear in the parameters windows once ypu run the
program.  Their only advantage is that in a problem with many defined
constants, of which only a few can be freely chosen, the parameter
window is not cluttered by dozens of parameters.  

Wiener parameters are more properly ``functions'' that return
scaled white noise. They are held fixed for $t$ to
$t+dt$ during an integration.  At each time step, they are then
changed and their value is a normally distributed random number with
zero mean and unit variance.  The program scales them by the
appropriate time step as well.  Their purpose is so that one can use
methods other than Euler for solving noisy problems.  In particular,
large steps can be taken using backward Euler without loss of
stability. 

The declarations:
\begin{verbatim}
par a[1..5]=.25
wiener w[1..5]
number z[1..5]=.123456
\end{verbatim}
behave in the obvious fashion.  
Note that this expression will lead to an error:
\begin{verbatim}
par a[1..2]=.5, c=.1234
\end{verbatim}
as XPP will expand it into 2 lines:
\begin{verbatim}
par a1=.5, c=.1234
par a2=.5, c=.1234
\end{verbatim}
which will give an ``duplicate name'' error.  On the other hand,
this expression will work:
\begin{verbatim}
par a[1..2]=.25,b[j]=.3
\end{verbatim}
and is the same as:
\begin{verbatim}
par a1=.25,b1=.3
par a2=.25,b2=.3
\end{verbatim}

 
\TCC{User-defined functions}
  User defined functions  have the following form:
\bvb
<name>(x1,x2,...)=<expression>
\end{verbatim}
where {\tt <name>} is the name of the function and {\tt x1,x2,...} are
the dummy arguments and {\tt <expression>} is a formula defining the
function. There can be at most 9 arguments.  



\tcc{Tables} Tables are another type of function but (at least as of
now) are only of one argument. 
The ``table'' declaration takes one of two forms: (i) file based and
(ii) function-based.  The file based version has the form:
\bvb
table <name> <filename>
\end{verbatim}
allows you to declare a function called {\tt
<name>} that reads in values from the file, {\tt <filename>} or as a
function of one variable over some interval and then
is used in your program as a function of 1 variable interpolated from
the tabulated values.  The values of this table are assumed to be
equally spaced and the file is an ASCII file with the format:
\bvb
<number of values>
<xlo>
<xhi>
y1
y2
.
.
.
yn
\end{verbatim}
Thus, $f(xlo)=y1$ and $f(xhi)=yn.$ (If the number of points in the
file description of the table starts 
with the ASCII character `i' , (e.g. i50 instead of 50)
then the interpolation will be piecewise
constant; otherwise it is linear.)
These tables can also be read in
from within XPP but a valid table must be given to start the program.
This table can be arbitrarily long (as memory permits) and thus you
can use experimental data as inputs to differential equations or even
sketch curves and use that as your nonlinearity.  

You can directly
input tabulated functions as well to speed up computations with
complicated functions.  In this case, after the table name put a
parenthesis symbol followed by the number of points, the minimum
argument and the maximum argument and then the function.  This should
be written as a function of ``t''.  Thus, the statement
\begin{verbatim}
table f % 501 -10 10 tanh(t)
\end{verbatim}
will produce a table of the hyperbolic tangent function from -10 to 10
consisting of 501 points.  

Using the data browser, you can create tabulated data from a
simulation to use later in a different simulation as a function or as
input or whatever.

\TCC{Global Flags}
The ``global'' declaration allows you to set some conditions and then
if these conditions hold reset dynamic variables according to the
conditions.  The form of ``global'' declarations is:
\bvb
global sign {condition} {name1=form1;...}
\end{verbatim}
The {\tt condition} is any combination of variables such
that if it is zero, then the desired event has occurred. 
{\tt sign} is either 1,-1, or 0.  A sign of 1 means that if the
condition goes from less than zero to greater than zero, the event has
occurred. A sign of -1 means that if the condition {\em decreases}
through zero that event has occurred.  Finally, a sign of 0, means
that any crossing of zero signals an event.  Each time the condition
is met, events occur.  There can be up to 20 events per condition and
they are delimited by braces and separated by semicolons.  Events are
always of the form: {\tt variable=formula} where {\tt variable} is one
of the differential equation variables and formula is some formula
involving the variables.  All formulae are first evaluated and then
the variables are updated.  Some examples are shown below.{\bf
WARNING!} The global flags are {\em ignored} by the adams
integrator.  Use GEAR, EULER, RUNGE-KUTTA, BACKWARD EULER, MODIFIED
EULER, STIFF, or QUALITY-RK.

{\em Note.} You will sometimes get the rather obscure message that the
program is ``Working too hard.''  This is a diagnostic that one of two
things is occurring.  First, due to round-off, sometimes the
extrapolation to zero for the condition is actually not zero but is
instead some very small number.  Thus XPP checks for the definition of
this which is called {\tt s } and if it is less than a user defined
value of {\tt smin} (changed in the numerics menu under ``sIng pt ctrl'')
then it is treated as zero.  The console reports this number when the
error message occurs so that you can change {\tt smin} to be larger,
say 1e-13, to avoid this message.  

The second situation that causes this to arise is that the time step
{\tt DT} is too large and the same condition is occuring twice in that
step.  Since the interpolation is linear, this means you should try to
take smaller steps; otherwise numerical errors will accumulate.  

Global flags can also be defined in groups, for example:
\begin{verbatim}
global 1 x[1..5]-1 {x[j]=0}
\end{verbatim}

\TCC{Initial and Boundary conditions}
The default for initial conditions of Markov and continuous variables
is zero.  Initial data can also be set in the ODE file (and of course
easily set from within XPP) in one of two ways:
\bvb
<name>(0)=value
\end{verbatim}
will set the variable {\tt <name>} to the specified
value. Alternatively, you can initialize many variables on one line by typing
\bvb
init <name1>=value1, <name2>=value2, ... 
\end{verbatim}

Boundary conditions can be placed anywhere in the file or ignored
altogether if you don't plan on solving boundary value problems. They
have the form:
\bvb
bndry <expression>
\end{verbatim}
where  stands for boundary condition. The expression is one
 involving your variables and which will be set to zero. In order to
 distinguish left and right boundary conditions, the following notation 
is used.  For the values of the variables at the left end of the
interval,
 use the symbol for
 the variable.  For the values at the right end of the interval, use the
 symbol for the variable appended by a single quote, '.  Thus to
 specify the boundary conditions $x(0)=1$, $y(1)=2$, the following is used:
\bvb
bndry x-1
bndry y'-2
\end{verbatim}
Periodic boundary conditions would be written as:
\bvb
bndry x-x'
bndry y-y'
\end{verbatim}
 Note that it is not necessary to specify the BCs at this point.  They can
 be specified within the program. Also note that the notation I have
 used allows the specification of mixed boundary conditions such as periodic
 BCs. (XPP has some additional special commands for periodic boundary
 conditions that enable the user to specify fewer equations than usual.)
 The number of boundary conditions must match the number of variables in
 your problem. Do not use the boundary value solver with Volterra,
delay, stochastic, or discrete equations.

Initial and boundary conditions can also be defined {\em en masse}
via:
\begin{verbatim}
bdry x[1..2]-2
x[1..2](0)=.345
\end{verbatim}
and this will be expanded in the expected fashion.

\TCC{Setting internal options} XPP has many many internal parameters
that you can set from within the program and four parameters that can
only be set before it is run.  Most of these internal parameters can be
set from the ``Options'' files described above and whose format is at
the end of this document.  However, it is often useful to put the
options right into the ODE file.  {\em NOTE: Any options defined in
the ODE file override all others such as the onres in the OPTIONS
file.} In addition, there are several options not available in the
options file.  These options are used by the ``silent'' integrator to
produce a file for output when running without X.

The format for changing the options is:
\bvb
@ name1=value1, name2=value2, ...
\end{verbatim}
where {\tt name} is one of the following and {\tt value} is either an
integer, floating point, or string.  (All names can be upper or lower
case). The first four options {\em can
only be set outside the program.}  They are:
\begin{itemize}
\item MAXSTOR={\tt integer} sets the total number of time steps that
will be kept in memory.  The default is 5000.  If you want to perform 
very long integrations change this to some large number.  
\item BACK= {\tt \{Black,White\}} sets the background to black or white.
\item SMALL={\tt fontname} where {\tt fontname} is some font available
to your X-server.  This sets the ``small'' font which is used in the
Data Browser and in some other windows.  
\item BIG={\tt fontname} sets the font for all the menus and popups.  
\end{itemize}

The remaining options can be set from within the program. They are 
\begin{itemize}
 \item XP=name sets the name of the variable to plot on the x-axis.
The default is {\tt T}, the time-variable.
\item YP=name sets the name of the variable on the y-axis.
\item ZP=name sets the name of the variable on the z-axis (if the plot
is 3D.) 
\item AXES={\tt \{2,3\}} determine whether a 2D or 3D plot will be
displayed.
\item TOTAL=value sets the total amount of time to integrate the
equations (default is 20).
\item DT=value sets the time step for the integrator (default is 0.05).
\item NJMP={\tt integer} tells XPP how frequently to output the
solution to the ODE.  The default is 1, which means at each
integration step.
\item T0=value sets the starting time (default is 0). 
\item TRANS=value tells XPP to integrate until {\tt T=TRANS} and then
start plotting solutions (default is 0.)
\item NMESH={\tt integer} sets the mesh size for computing nullclines
(default is 40).
\item METH={\tt \{
discrete,euler,modeuler,rungekutta,adams,gear,volterra, backeul,
qualrk,stiff,cvode\}}
sets the integration method (see below; default is Runge-Kutta.)
\item DTMIN=value sets the minimum allowable timestep for the Gear
integrator.
\item DTMAX=value sets the maximum allowable timestep for the Gear
integrator
\item ATOLER=value sets the absolute tolerance for CVODE.
\item TOLER=value sets the error tolerance for the Gear, adaptive RK,
and stiff integrators. It is the relative tolerance for CVODE. 
\item BOUND=value sets the maximum bound any plotted variable can
reach in magnitude. If any plottable quantity exceeds this, the
integrator will halt with a warning.  The program will not stop
however (default is 100.)
\item DELAY=value sets the maximum delay allowed in the integration
(default is 0.)
\item XLO=value,YLO=value,XHI=value,YHI=value set the limits for
two-dimensional plots (defaults are 0,-2,20,2 respectively.) Note that
for three-dimensional plots, the plot is scaled to a cube with
vertices that are $\pm1$ and this cube is rotated and projected onto
the plane so setting these to $\pm2$ works well for 3D plots.
\item
XMAX=value, XMIN=value, YMAX=value, YMIN=value, ZMAX=value, ZMIN=value set
the scaling for three-d plots.
\item OUTPUT=filename sets the filename to which you want to write for
``silent'' integration.  The default is ``output.dat''. 
\item POIMAP={\tt \{ section,maxmin\} } sets up a Poincare map for
either sections of a variable or the extrema.  
\item POIVAR=name sets the variable name whose section you are
interested in finding.
\item POIPLN=value is the value of the section; it is a floating
point.
\item POISGN={\tt \{ 1, -1, 0 \}} determines the direction of the
section.  
\item POISTOP=1 means to stop the integration when the section is
reached.
\item RANGE=1 means that you want to run a range integration (in batch
mode). 
\item RANGEOVER=name, RANGESTEP, RANGELOW, RANGEHIGH, RANGERESET={\tt
Yes,No}, RANGEOLDIC={\tt Yes,No} all correspond to the entries in the
range integration option (see below).
\item TOR\_PER=value, defined the period for a toroidal phasespace and
tellx XPP that there will be some variables on the circle.
\item FOLD=name, tells XPP that the variable <name> is to be
considered modulo the period.  You can repeat this for many variables.
\end{itemize}

\TCC{Named parameter sets}  Sometimes, you wnat to prepare a bunch
of simulations that use different initial data or parameter values or
numerical methods, etc.  You can, of course, change these within the
program by choosing the desired option and changing it.  Or, you can
do the simulation and save the results in a ``.set'' file (see
below). Another way to do this is by adding a bunch of parameter sets
to the ode file. The format for this is:
\begin{verbatim}
set name {item1=value1, item2=value2, ..., }
\end{verbatim}
Then, you tell XPP to use the named set and it will do all of the
things inside the brackets.  The items are any parameter name, any
variable name, or any of the internal options named above. The named
sets are accessed through the {\tt Get par set} menu item in the {\tt
File} menu.
Here is an example:
\begin{verbatim}
# test
x'=-a*x+c
par a=1,c=1
set set1 {a=1,c=1,x=0,dt=.25}
set set2 {a=.25,c=0,x=1,dt=.1}
done
\end{verbatim}
If you load ``set1'' then the parameters, initial conditions, and
{\tt DeltaT} will be set to the values in the brackets.  Choosing
``set2'' sets them differently.  The names can be anything you like.

\TCC{Finishing up}
The last line in the file should be ``done'' telling the ODE reader
that the file is over.


\TCC{Important shortcut} All of the declarations, {\tt
markov, parameter, wiener, table, aux, init, bndry, global, done, }
can be abbreviated by their first letter.

\tc{Reserved words} 
\begin{center}
\Large Reserved words 
\end{center}
You should be aware of the following
keywords that should not be used in your ODE files for anything other
than their meaning here.
\bvb
sin cos tan atan atan2 sinh cosh tanh
exp delay ln log log10 t pi if then else
asin acos heav sign ceil flr ran abs 
max min normal besselj bessely erf erfc
arg1 ... arg9  @ $ + - / * ^ ** shift
| > < == >= <= != not \# int sum of i'
\end{verbatim}




These are mainly self-explanatory. The nonobvious ones are:
\begin{itemize}
\item {\tt heav(arg1)} the step function, zero if {\tt arg1<0} and 1 otherwise.
\item {\tt sign(arg)} which is the sign of the argument (zero has sign 0)
\item {\tt ran(arg)} produces a uniformly distributed random number
between 0 and {\tt arg.}
\item {\tt besselj, bessely } take two arguments, $n,x$ and return
respectively, $J_n(x)$ and $Y_n(x),$ the Bessel functions.
\item { \tt erf(x), erfc(x)} are the error function and the
complementary function. 
\item {\tt normal(arg1,arg2)} produces a normally distributed random number
with mean {\tt arg1}  and variance {\tt arg2}.
\item {\tt max(arg1,arg2)} produces the maximum of the two arguments
and {\tt min}  is 
the minimum of them.
\item {\tt if(<exp1>)then(<exp2>)else(<exp3>)} evaluates {\tt <exp1> }
If it is nonzero 
it evaluates to {\tt <exp2>} otherwise it is { \tt <exp3>}.  E.g. {\tt if(x>1)then(ln(x))else(x-1)}
will lead to {\tt ln(2)}  if {\tt x=2}  and { \tt -1 if x=0.}
\item {\tt delay(<var>,<exp>)} returns variable {\tt <var>} delayed by the result of
 evaluating {\tt <exp>}.  In order to use the delay you must inform
the program of the maximal possible delay so it can allocate storage.
(See the section on the NUMERICS menus.)
\item {\tt  ceil(arg),flr(arg)}  are the integer parts of{\tt  <arg>} returning the
 smallest integer greater than and the largest integer less than {\tt <arg>}.  
\item  {\tt t } is the current time in the integration of the differential equation.
\item {\tt  pi}  is $\pi.$ 
\item {\tt arg1, ..., arg9} are the formal arguments for functions 
\item {\tt int, \#} concern Volterra equations.
\item {\tt shift(<var>,<exp>)} This operator evaluates the expression
{\tt <exp>} converts it to an integer and then uses this to indirectly
address a variable whose address is that of {\tt <var>} plus the
integer value of the expression.  This is a way to imitate arrays in
XPP.  For example if you defined the sequence of 5 variables, {\tt
u0,u1,u2,u3,u4} one right after another, then {\tt shift(u0,2)} would
return the value of {\tt u2.} 
\item {\tt sum(<ex1>,<ex2>)of(<ex3>)} is a way of summing up things.
The expressions {\tt <ex1>,<ex1>} are evaluated and their integer
parts are used as the lower and upper limits of the sum.  The index of
the sum is {\tt i'} so that you cannot have double sums since there is
only one index.  {\tt <ex3>} is the expression to be summed and will
generally involve {\tt i'.}  For example  {\tt sum(1,10)of(i')} will
be evaluated to 55.  Another example combines the sum with the shift
operator.  {\tt sum(0,4)of(shift(u0,i'))} will sum up {\tt u0} and the
next four variables that were defined after it.  An example below
shows how this can be used to solve large systems of equations that
are densely coupled. 
\end{itemize}
\TCC{Compatibility} The parser in this distribution is a great
improvement over the old style parser. 
The parser distinguishes between the old style and the
new style by whether or not the first line of the file contains a
number.  If the first line is a number, then the old style parser is
used.  Otherwise, the new style is used.  
\section{Command line arguments}
There are only (at present) three command line arguments.
They are 
\begin{description}
\item{-xorfix} This changes the way rubber-band drawing for zooms and
other things is done. If you do not see a box when you zoom in, you
should run XPP with this argument.
\item{-convert}  This allows you to convert old style parser format to
the new style which is much more readable. The program creates a file
with the same name as the input file but with the extension {\tt .new}
appended. It works on on the examples I have tried but it is still in
beta testing.
\item{-silent} This allows you to run XPP's integrators without using
the X-windows stuff.  The result of the integration is saved to a file
called ``output.dat'' but this can be changed. The length of
integration, methods, Poincare sections, etc, are all specified in
either the options file (see section \ref{optfiles})
 or in the internal options.
\end{description}
The only other thing on the command line should be the file name.  
Thus, 
\begin{verbatim}
xppaut test.ode -xorfix -convert
\end{verbatim}
will convert {\tt bob.ode} to the new format and run it with the {\tt
xorfix}.

\section{ Examples}

Nothing helps one understand how to use a program better than lots of
examples.  
 
\tc{Morris-Lecar}
\begin{center}
Morris-Lecar Equations
\end{center}
The Morris-Lecar equations arise as a simplification of a model for
barnacle muscle oscillations.  They have the form:
\beqann
C\frac{dV}{dt} &=&
g_L(V_L-V)+g_Kw(V_k-V)+g_{Ca}m_\infty(V)(V_{Ca}-V)+I \\
\frac{dw}{dt} &=& \phi \lambda_w(V)(w_\infty(V)-w)
\eeqann
where
\beqann
m_\infty(V) &=& .5(1+\tanh((V-V_1)/V_2)) \\
w_\infty(V) &=& .5(1+\tanh((V-V_3)/V_4)) \\
\lambda_w &=& \cosh((V-V3)/(2V_4))
\eeqann
This is a very straightforward model and so the equation file is
pretty simple:
\begin{verbatim}
# The Morris-Lecar equations 

# Declare the parameters
p gl=.5,gca=1,gk=2
p vk=-.7,vl=-.5,vca=1
p v1=.01,v2=.145,v3=.1,v4=.15
p i=.2,phi=.333   
  
# Define some functions
minf(v)=.5*(1+tanh((v-v1)/v2))
winf(v)= .5*(1+tanh((v-v3)/v4))
lamw(v)= cosh((v-v3)/(2*v4))

# define the right-hand sides
v'= gl*(vl-v)+gk*w*(vk-v)+gca*minf(v)*(vca-v)+i
w'= phi*lamw(v)*(winf(V)-w)

# some initial conditions -- not necessary but for completeness
v(0)=.05
w(0)=0

# Done!!
d
\end{verbatim}
Note that some errors are now caught by the parser.  For
example, duplicate names and illegal syntax are found.  


Suppose you want to keep track of the calcium current as an auxiliary
variable.  Then, the following file will work
\begin{verbatim}
# The Morris-Lecar equations

# Declare the parameters
p gl=.5,gca=1,gk=2
p vk=-.7,vl=-.5,vca=1
p v1=.01,v2=.145,v3=.1,v4=.15
p i=.2,phi=.333   
  
# Define some functions
minf(v)=.5*(1+tanh((v-v1)/v2))
winf(v)= .5*(1+tanh((v-v3)/v4))
lamw(v)= cosh((v-v3)/(2*v4))

# define the right-hand sides
v'= gl*(vl-v)+gk*w*(vk-v)-gca*minf(v)*(v-vca)+i
w'= phi*lamw(v)*(winf(v)-w)
#
aux ica=gca*minf(v)*(v-vca)

# some initial conditions -- not necessary but for completeness
v(0)=.05
w(0)=0

# Done!!
d
\end{verbatim}
Note that we are wasting computational time
since we compute $I_{Ca}$ twice; once as a contribution to the
potential change and once as an auxiliary variable.  In a FORTRAN or
C program, one would compute it as a local variable and use it in both
instances.  This is the purpose of fixed variables.  (For the present
problem, the computational overhead is trivial, but for coupled arrays
and other things, this can be quite substantial.)  The last program
uses a fixed variable to reduce the computation:
\begin{verbatim}
# The Morris-Lecar equations ml1.ode

# Declare the parameters
p gl=.5,gca=1,gk=2
p vk=-.7,vl=-.5,vca=1
p v1=.01,v2=.145,v3=.1,v4=.15
p i=.2,phi=.333   
  
# Define some functions
minf(v)=.5*(1+tanh((v-v1)/v2))
winf(v)= .5*(1+tanh((v-v3)/v4))
lamw(v)= cosh((v-v3)/(2*v4))

# define the right-hand sides
v'= gl*(vl-v)+gk*w*(vk-v)-icaf+i
w'= phi*lamw(v)*(winf(v)-w)

# where
icaf=gca*minf(v)*(v-vca)

# and
aux ica=icaf

# some initial conditions -- not necessary but for completeness
v(0)=.05
w(0)=0


# Done!!
d
\end{verbatim}
The calcium current is computed once instead of twice.  This is why
``fixed'' variables are useful.


{\bf NOTE:}  Since fixed quantities are not ``visible'' to the user
and auxiliary quantities are not ``visible'' to the internal formula
compiler, we can use the same name for the fixed as the
auxiliary variable; the {\tt aux } declaration essentially makes it
visible to the user with almost no computational overhead. However, it
does generate an error message (not fatal) so it is best to make all
names unique.  
   
\tc{Cable model} \begin{center}
A linear cable equation with boundary conditions
\end{center}
In studying a dendrite, one is often interested in the steady-state
voltage distribution.  Consider a case where the dendrite is held at
$V=V_0$ at $x=0$ and has a leaky boundary condition at $x=1.$
Then the equations are:
\beqann
\lambda^2\frac{d^2V}{dx^2} &=& V(x) \\
V(0)&=& V0 \\
a\frac{dV(1)}{dx} + b V(1) &=& 0
\eeqann
When $a=0,b\ne0$ the voltage at $x=1$ is held at 0.  When $a\ne0,b=0$
there is no leak from the cable and the conditions are for sealed end.
Since this is a second order equation and XPP can only handle first
order, we write it as a system of two first order equations in the
file which is:
\begin{verbatim}
# Linear Cable Model cable.ode

# define the 4 parameters, a,b,v0,lambda^2
p a=1,b=0,v0=1,lam2=1
# now do the right-hand sides
v'=vx
vx'=v/lam2

# The initial data
v(0)=1
vx(0)=0

# and finally, boundary conditions
# First we want V(0)-V0=0
b v-v0
#
# We also want aV(1)+bVX(1)=0
b a*v'+b*vx'
# Note that the primes tell XPP to evaluate at the right endpoint
d
\end{verbatim}
Note that I have initialized $V$ to be 1 which is the appropriate
value to agree with the first boundary condition.  

\tc{Delayed feedback}\begin{center}
The delayed inhibitory feedback net
\end{center}
A simple way to get oscillatory behavior is to introduce delayed
inhibition into a neural network.  The equation is:
\[
\frac{dx(t)}{dt} = -x(t) + f(ax(t)-bx(t-\tau)+P)
\]
where $f(u)=1/(1+\exp(-u))$ and $a,b,\tau$ are nonnegative parameters.
Here $p$ is an input.  The XPP file is:
\begin{verbatim}
# delayed feedback

# declare all the parameters, initializing the delay to 3
p tau=3,b=4.8,a=4,p=-.8
# define the nonlinearity
f(x)=1/(1+exp(-x))
# define the right-hand sides; delaying x by tau
dx/dt = -x + f(a*x-b*delay(x,tau)+p)
x(0)=1
# done
d
\end{verbatim}
Try this example after first going into the numerics menu and changing
the maximal delay from 0 to say, 10.  Glass and Mackey describe a few
other delay systems some of which have extremely complex behavior.
Try them out.



\tc{Random mutation}\begin{center}
A population problem with random mutation
\end{center}
This example illustrates the use of Markov variables and is due to Tom
Kepler.  There are two variables, $x_1,x_2$ and a Markov state
variable, $z.$  The equations are:
\beqann
x_1' &=& x_1(1-x_1-x_2) \\
x_2' &=& z(ax_2(1-x_1-x_2)+\epsilon x_1)
\eeqann
and $z$ switches from 0 to 1 proportionally to $x_1.$  The transition
matrix is 0 everywhere except in the 0 to 1 switch where it is
$\epsilon x_1.$ So z=1 is an absorbing state.

Initial conditions should be $x_1(0)=1.e-4$ or so, $x_2(0)$ the same (this is just
for convenience; we really want $x_2(0)=0$ and then have it jump discontinuously
to $1.e-4$ or so when $z$ makes its transition, but this shouldn't matter that much)


This models the population dynamics of two populations $x_1,x_2$ in competition
with each other.  $x_2$, initially absent, is a mutant of $x_1.$  The mutation rate
$\epsilon$ should be smaller than one.  The relative advantage, $a$,
should be larger than one.

One expects that $x_1$ grows for a while, eventually $z$ makes its transition,
$x_2$ begins to grow and eventually overtakes $x_1.$ 

The equation file for this is
\begin{verbatim}
# Kepler model kepler.ode

init x1=1.e-4,x2=1.e-4,z=0
p eps=.1,a=1
x1' = x1*(1-x1-x2)
x2'= z*(a*x2*(1-x1-x2)+eps*x1)
# the markov variable and its transition matrix
markov z 2
{0} {eps*x1}
{0} {0}
d
\end{verbatim}
\medskip

\tc{Flags and discontinuities}\begin{center}
Some equations with flags
\end{center}
Tyson describes a model which involves cell growth:
\beqann
\frac{du}{dt} &=& k_4(v-u)(a+u^2)-k_6u \\
\frac{dv}{dt} &=& k_1m-k_6u \\
\frac{dm}{dt} &=& bm 
\eeqann

This is a normal looking model with exponential growth of the mass.
However, if $u$ decreases through 0.2, then the ``cell'' divides in
half.  That is the mass is set to half of its value.  Thus, we want to
flag the event $u=.2.$  The file in the {\em new format} (no sense
using the old format since global variables did not appear in earlier
versions of XPP)
\begin{verbatim}
# tyson.ode 
i u=.0075,v=.48,m=1
p k1=.015,k4=200,k6=2,a=.0001,b=.005
u'=  k4*(v-u)*(a+u^2) - k6*u
v'= k1*m - k6*u
m'= b*m
global -1 {u-.2} {m=.5*m} 
d
\end{verbatim}
Everything is fairly straightforward.  When $u-.2$ decreases through
zero tha is, $u$ is greater than .2 and then less than .2, the mass is
cut in half.  Integration using any of the integrators {\em except
ADAMS!} yields a regular limit cycle oscillation.

Another example with discontinuities arises in the study of coupled
oscillators.  Two phase oscillators $x,y$ travel uniformly around the
circle.  If $x$ hits $2\pi$ it is  reset to 0 and adds an amount
$r(y)$ to the phase of $y.$  The XPP file is
\begin{verbatim}
# delta coupled oscillators delta.ode
x'=1
y'=w
global 1 {x-2*pi} {x=0;y=b*r(y)+y}
global 1 {y-2*pi} {y=0;x=b*r(x)+x}
r(x)=sin(x+phi)-sin(phi)
par b=-.25,phi=0,w=1.0
done
\end{verbatim}
Note that there are 2 conditions and each creates 2 events.
\medskip

\tc{Large coupled systems}\begin{center}
Large coupled systems
\end{center}
Many times you want to solve large coupled systems of differential
equations.  Writing the ODE files for these can be a tedious exercise.
The array declaration described above simplifies the creation of ODE
file for this.
For
example, suppose that you want to solve:
\[
\frac{du_j}{dt} = -u_j + f(a\sum_{i=0}^{n-1}\cos\beta(i-j) u_i)
\]
which is a discrete convolution.  Suppose that {\tt n=20} so that this
is 20 equations.  Its not very convenient to type these equations so
instead you can create a simple file in which the equations are typed
just once:
\begin{verbatim}
# chain of 20 neurons
param a=.25,beta=.31415926
u[0..19]'=-u[j]+f(a*sum(0,19)of(cos(beta*([j]-i'))*shift(u0,i')))
f(x)=tanh(x)
done
\end{verbatim}
The {\tt u[0..19]} line with an index spanning {\tt 0} to {\tt 19}
tells XPP to repeat this from 0 to 19.  Thereafter, the string {\tt
[j]} or some simple arithmetic expressions of {\tt j} are evaluated
and substituted verbatim.  
Note how I have used the shift operator on {\tt u0} to make it act like
an array. (Recall that {\tt shift(x,n)} gives the value of the
variable that is defined {\tt n} after the variable {\tt x} is
defined.  

Here is one more example of nearest neighbor coupling in an excitable
medium.  This is a discretization of a PDE.  I will use the analogue
of ``no flux'' boundaries so that the end ODEs will be defined separately.  
\begin{verbatim}
# pulse wave
param a=.1,d=.2,eps=.05,gamma=0,i=0
v0(0)=1
v1(0)=1
f(v)=-v+heav(v-a)
#
v0'=i+f(v0)-w0+d*(v1-v0)
v[1..19]'=i+f(v[j])-w[j]+d*(v[j-1]-2*v[j]+v[j+1])
v20'=i+f(v20)-w20+d*(v19-v20)
#
w[0..20]'=eps*(v[j]-gamma*w[j])
@ meth=modeuler,dt=.1,total=100,xhi=100
done
\end{verbatim}

I only need to define {\tt v0,v20} separately since {\tt w} does not diffuse.
Note how much typing I saved.  




\bigskip











 
  
XPP comes with some other examples that I urge you to look at. A
Volterra example is shown below.

\section{Using xpp}
XPP runs from the command line by typing
\bvb
 xpp <filename>
\end{verbatim}
where {\tt <filename>} is an ODE file as described above.  This argument is optional. 
 If you invoke XPP without the argument, you are asked whether you want to
 read or create a file. Choosing {\tt (r)} to read the file prompts you for an ODE
 file name.  You can get a directory by pressing {\tt <RETURN>}.  If you choose to 
create a new file, you must write a set of commands such as above and then
 save the file.  	It will bounce you out if you make any mistakes. 
 It is not the recommended way to do it -- use an editor and read in the file. 
 XPP will look at the extension of a file and if it is {\tt .dif } the program
 will assume that the model is a difference equation and will treat it as a map. 
 The solution, equilibria, stability and listing of difference equations 
is treated differently from differential equations.  The default plot style 
is also plotting of points instead of lines as for the continuous time case.

\tc{Windows}
\tcc{Main Window}

Once you get a file loaded and it is valid, the program and all of the 
windows will pop up.  There is the main window with the graphics and menus.  
This has several regions on it: the graphics window, the menu window, the command
 window across the top, the message window at the bottom, and three
sliders for binding to parameters.  Initial 
conditions and parameters as well as text labels can be assigned
values in the command window for those who don't like the mouse.  The message 
window gives the position of the pointer and ``tips'' about the
various menu items.  If 
you click the mouse in the main graph and hold it down as you move it, the
coordinates appear in the message window.  The left 
most region is the menu.  Pressing the mouse on a command is the same 
as typing the hot letter for the command and will execute it.  As you
move through the menu items, a brief hint about its maning appears in
the message window. The pointer is 
a hand in the menu window.  {\em In addition to the main window, there
several other windows that will be created.  Due to user demand, these
are all instantly iconified so that you must open them up if you want
to use them.}

\tcc{Initial data,delays,parameters}

In addition to the main window, there is the initial data window which
 contains the names of the variables and the current initial data.  The 
delay window is similar containing the initial data for delay equations. 
 The parameter window contains the list of all parameters and their current
 values.  The boundary condition window is similar with the list of boundary
 conditions given.  These 4 windows constitute the main user editable windows.
  All of the quantities and formulae contained within these can be edited within
 the program.  Clicking on any item in these windows will cause the
 cursor to appear where it can be edited.  If you type the {\tt Enter}
 key after entering the item, the cursor moves to the next item and
 the edited item is accepted.  You can press the {\tt OK} button also
 to accept the item and the pointer is then redirected to the main
 window.  Pressing the {\tt Default} key in the Parameter
 Window will reset the parameters to their values when the program
 boots up.  The {\tt Cancel} button doen't do much of anything yet,
 but maybe someday if I can think of a use! 
 Pressing the {\tt Tab} key is the same as pressing the {\tt Ok} button.
 Pressing {\tt Esc}  exits without changing. 

You can also press the letter {\tt p}  in the main window
to change parameters or simply click
 on the parameter command in the main menu.  Typing {\tt default } will load all 
the parameters with their values when the program first fires up. The
initial data can also be changed from the command window (see below.)  

\tcc{Equation window}

Another window contains the listing of the differential equations. If
 the file is a difference equation, the listing shows this.  These equations 
can be scrolled up and down by pressing the mouse on the {\tt <Up> <Down>} buttons
 or using the cursor keys. Resizing lets you see more and closing iconifies the
 window. Once you are in XPP you cannot alter the differential equations.  Thus
, to solve a different problem you must exit and start again.  I hope to change 
this in the future.  All of these windows can be iconified with no effect on the
 performance of the program.   

\tcc{Browse window}

The other major window contains listings of the numerical values of the solutions 
to the differential equations. 
 This browser window allows you to postprocess
 data and save and load numbers to the disk.
The first few lines contain a variety
of buttons.  The next line has hints about the various buttons that
are given as you move the mouse over the buttons. The names of the
variables, time, and auxilliary variabls are next.  The the numbers
themselves are given. 
  Pressing the cursor keys, the
 {\tt HOME, END, PgUp,PgDn}  keys scrolls you through tthe data either up and down or
 left and right.  Resizing the window lets you see more data.  There are many
 other keys and you can click on the buttons or use the keyboard if the mouse
 pointer is in the Browser.  You can load or save data files by choosing the 
{\tt Load}  and {\tt Save}  keys.  You will be asked for a filename which you should enter.
  Choosing {\tt Get} will load the top row of numbers into the initial data.
 {\tt Find} will find the data row for which the specified variable is closest to
 the specified value.   {\tt Replace} allows you to replace a column by some formula 
involving the column itself and all other columns.  You are asked for the column 
to replace.  Type in one of the names across the top row.  Then you are asked 
for the formula.  Type it in and the entire column will be replaced.  For example 
if the names of your variables are {\tt X}  and { \tt Y}  and you
choose {\tt X } as the column
 to replace and type 
\bvb 
       X^2 + T*Y
\end{verbatim}
then the column {\tt X } will be replaced by the function at each row.  If the first 
character in the formula is the {\tt \& } symbol, an integral of the formula will be
 placed in the specified column.  If the formula is anything of the form :
\bvb
	@<variable>
\end{verbatim}
the numerical derivative of the {\tt <variable> } will be computed ad placed in the
 column. 
Choosing {\tt Unreplace} will restore the most recent column.  Each time you integrate
 the equations or solve a boundary value problem or look at a range of equilibria, 
the numbers are replaced by the most recent.  You must save the data to a
 file to use it again.  Pressing {\tt First}  and {\tt last} will mark the data so that you 
can save a portion or restore a portion.  {\tt Home} takes you to the beginning and
 {\tt End} to the end of the list. 

The {\tt Table} key allows you to save a column as a ``table'' file
which can be used in XPP as a function.  

\tcc{Parameter sliders}  These are really cool and I stole the idea
from Madonna, a Mac ODE solver. Click on one of them where it says
{\tt Parameter} and fill it in.  These bind the slider to a particular
parameter.  You set the low and high values of the slider and the
current value.  Click {\tt OK} to accept. Change them by clicking on
the name window.  Then move the mouse to the slider and pressing it,
move it around.  When done, you have set the parameter to the
indicated value.


Other windows are opened on occasion and can remain open throughout the
 session.  The usual methods for resizing windows works on the graphics and
 listing and browser windows.

When there are multiple graphics windows, only one is active at a given time
 and only one can be drawn to.  The window is made active by clicking the
 left button when inside the window.  A little square appears in the upper left
 corner of the window to indicate its activity.  The variables on the
 axes are given in the title of the window and are Z-Axis vs Y-Axis vs X-Axis.  


Two other windows appear permanently, when needed: the Equilibrium window
 that gives values and stability of the most recently computed
 equilibrium and the AUTO window for bifurcation analysis.


The main commands are on the first menu and the numerics menu.  Communication 
to the program is via the keyboard and the mouse and through parameter boxes
 that are displayed.  Almost all commands have keyboard shortcuts so that if
 you know the original DOS version, you can use the X version with little help.  
When you are prompted for floating point numbers at the command line, then 
you can insert either numbers or mathematical expressions provided the first
 character is the {\tt \%} symbol. Thus if you want to input {\tt sin(1.5)} as a floating 
point number, you would write
\bvb
          %sin(1.5)
\end{verbatim}
at the command line. You can enter such expressions in the parameter and
initial data windows and these are converted to numbers.  
 

\section{The main commands}
All commands can be invoked by typing the hot key for that command (capitalized
 on the menu) or clicking on the menu with the mouse.  Usually most commands can
 be aborted by pressing the {\tt Esc} key. Once one of these is chosen, the
 program begins to calculate and draw the trajectories.  If you want to
 stop prematurely, press the {\tt Esc} key and the integration will
stop. 

\begin{description}

\tc{ICs}\item[(I)nitial conds] This invokes a list of options for integrating the
 differential equations.  The choices are:
\begin{description}
\item[(G)o] which uses the initial data in the IC window and the current
 numerics parameters to solve the equation.  The output is drawn in the
 current selected graphics window and the data are saved for later use.  The
 solution continues until either the user aborts by pressing {\tt Esc}, the
 integration is complete, or storage runs out.
\item[(N)ew] This prompts you at the command line for each initial condition.
 Press {\tt Return} to accept the value presented.
\item[(L)ast]  This uses the end result of the most recent integration as the
 starting point of the curret integration.
\item[(S)hift] This is like Last except that the stating time is shifted to the 
current integration time.  This is irrelevant for autonomous systems but is
 useful for nonautonomous ODEs.
\item[(O)ld] This uses the most recent initial data as the current initial 
data.  It is essentially the same as Go.
\item[(R)ange]  This lets you integrate multiple times with the results shown
 in the graphics window.  Pressing this option produces a new window with
 several boxes to fill in.  First choose the quantity you want to range
 over.  It can be a parameter or a variable.  The integrator will be called
 and this quantity will be changed at the beginning of each integration.
  Then choose the starting and ending value and the number of steps. The
 option Reset storage only stores the last integration.  If you choose not
 to reset, each integration is appended to storage.  Most likely, storage
 will be exceeded and the integration will overwrite or stop. The option 
to use last initial conditions will automatically use the final result of
 the previous integration as initial dat for the next integration.  Otherwise,
 the current ICs will be used at each step (except of course for the  variable
 through which you are ranging.) When you are happy with the parameters, simply
 press the OK button.  Otherwise, press the Cancel button to abort.  Assuming
 that you have accepted, the program will compute the trajectories and plot 
them storing none of them or all of them.  If you press {\tt Esc} it will abort the 
current trajectory and move on to the next.  Pressing the {\tt / } key will abort
 the whole process. If you choose {\tt Yes} in the {\tt Movie} item,
then after each integration, XPP will take a snapshot of the picture.
You can then replay this series of snapshots back using the Kinescope. 
\item[(M)ouse]  allows you to specify the values with the mouse.  Click at the
 desired spot.
\item[s(H)oot] allows you to use initial data that was produced when
you last searched for an equilibrium.  When a rest state has a single
positive or negative eigenvalue, then XPP will ask if you want to
approximate the invariant manifold.  If you choose {\tt yes} to this,
then the initial data that were used to compute the trajectories are
remembered.  Thus, when you choose this option, you will be asked for
a number 1-4.  This number is the order in which the invariant
trajectories were computed.  Note if the invariant set is a stable
manifold, then you should integrate backwards in time.
\item[(F)unction] allows you to set all the initial data as a
formula. This is good for systems that represent chains of many ODEs.
The formula is given as a function of {\em t } and the initial data of
the $i^{th}$ variable is given by evaluating this formula at $t=i.$
Thus, {\tt 6.28*ran(1) } will randomly set each variable to some
number between 0 and 6.28.  You will be asked for the name of the
first variable for the formula, then number of variables to set and
the formula.   If you just type a carriage return for the variable
name, the program will stop asking. Thus for example, if you have a
system with two arrays {\tt u0,...,u10} and {\tt v0,...,v10}, then
when asked for the first element, type {\tt u0} and then {\tt 11} for
the number of elements, then type in a formula of {\tt t} for the {\tt
thh} element in the array.  Then when prompted again, type in {\tt v0}
etc. When prompted a third time, type {\tt Enter} and the integration
will begin again.  
  
\end{description}
\tc{Continue}\item[(C)ontinue]  This allows you to continue integrating appending the data
 to the 
current curve. Type in the new ending time.
\tc{Nullclines}\item[(N)ullclines]  This option allows you to draw the nullclines of systems.  
They are most useful for two-dimensional models, but XPP lets you draw them for
 any model.  The constraints are the same as in the direction fields option
 above.  The menu has 4 items.
\begin{description}
\item[(N)ew] draws a new set of nullclines.
\item[(R)estore] restores the most recently computed set.
\item[(A)uto] turns on a flag that makes 
XPP redraw them every time it is necessary because some other window obscured
 them.
\item[(M)anual] turns this flag off so that you must restore them manually.
 The X-axis nullcline is blue and the Y-axis nullcline is red.
\end{description}
\tc{Direction Field/Flow}\item[(D)irection Field/Flow]
This option is best used for two-dimensional systems however,
 it can be applied to any system.  The current graphics view must be a
 two-d plot in which both variables are different and neither is the time
 variable, T.  There are two items.
\begin{description}
\item[(D)irection fields] Choosing the direction field option will prompt you for 
a grid size.  The two dimensional plane is broken into a grid of the size
 specified and lines are drwan at each point specifying the direction of 
the flow at that point.  If the system is more than two-dimensional, the
 other variables will be held at the values in the initial conditions
window.
\item[(F)low] Choosing the flow, you will be prompted for a grid size
and
 trajectories
 started at each point on the grid will be integrated according to the
 numerical parameters.  Any given trajectory can be aborted by pressing 
{\tt Esc} and the whole process stopped by pressing {\tt /}.  The
remaining
 variables 
if in more than two-dimensions are initialized with the values in the
IC window. 
\item[(N)o Dir Field] turs off redrawing of direction fields when you
click on (Redraw).  Erasing the screen automatically turns this off.
\item[(C)olorize] This draws a grid of filled rectangles on the screen
whose color is coded by the velocity or some other quantity. (See the
Numerics Colorize menu item).
\end{description}


\tc{Window}\item[(W)indow] allows you to rewindow the current graph.  Pressing this presents 
another menu with the choices:
\begin{description}
	\item[(W)indow]  A parameter box pops up prompting you for the values.  
Press OK or CANCEL when done.
	\item[(Z)oom in]  Use the mouse to expand a region by clicking, dragging and 
releasing.  The view in the rectangle will be expanded to the whole window.
	\item[Zoom (O)ut]  As above but the whole window will be shrunk into the 
rectangle.
	\item[(F)it]  The most common command will automatically fit the window so 
the entire curve is contained within it.  For three-D stuff the window data 
will be scaled to fit into a cube and the cube scaled to fit in the window. 
Use this often.
\end{description}
(NOTE:  On some displays, no rubber boxes are drawn for the zooming
operations.  If this occurs, define the {\tt XORFLAG} and recompile
the program.  See the Makefile.)
\tc{Phase Space}\item[ph(A)se space]  XPP allows for periodic domains so that you can solve 
equations on a torus or cylinder.  You will be prompted to make (A)ll variables
periodic, (N)o variables periodic or (C)hoose the ones you want.  You will 
be asked for the period which is the same for all periodic variables (if they 
must be different, rescale them) Choose them by name or by clicking the 
appropriate names from the list presented to you. An {\tt X} will
appear next to the selected ones.  Clicking toggles the {\tt X}. 
 Type {\tt Esc} when done or CANCEL or DONE.
 XPP mods your variables by this period and is smart
 enough when plotting to not join the two ends. 

\tc{Kinescope}\item[(K)inescope]  This allows you to capture the bitmap of the active window and
 play it back.  It works poorly in X, you cannot save the bitmaps, nor will
 it run properly across the network.  This is not my fault, it is idiosyncratic
 with X.  It is slow also. On my HP, however, it is quite smooth so I
guess it depends on the graphics processors. 
  My PC makes much smoother animations. In any case, 
another menu pops up with the choices:
\begin{description}
\item[(C)apture] which takes a snapshot of the 
currently active window
\item[(K)ill] which deletes all the snapshots
\item[(A)utoplay] which 
continuously plays back snapshots until you press the middle mouse
button
\item[(P)layback] which cycles thru the pictures each time you click the left mouse 
button and stops if you click the middle.  It will beep when done.
\end{description}


\tc{Graphics}\item[(G)raphic stuff]  This induces a popup menu with several
choices.
\begin{description}
\item[(A)dd curve]  This lets you add another curve to the picture. A 
parameter box will appear aking you for the variables on each axis, a color,
 and line type (only 1 - solid and 0 - dotted are allowed currently) All 
subsequent integrations and restorations will include the new graph.  Up to 
10 per window are allowed.
\item[(D)elete last]  Will remove most recent curve.
\item[(A)ll but first] Deletes all curves but first.
\item[(E)dit curve]  You will be asked for th curve to edit.  The first is 0,
 the second 1, etc.  You will get a parameter box like the add curve
option.
\item[(P)ostscript]  This will ask you for a file name and write a postscript
 representation of the current window.  Nullclines, text, and all
graphs will be 
 plotted.  You will be asked for Black and White or Color.  Color tries
to match the color on the screen (but will not draw color-coded
trajectories yet). Black and white will use a variety of dashed curves
to plot the curves.  
\item[(F)reeze]  This will create a permanent curve in the window.  Usually, 
when you reintegrate the equations or load in some new data, the current curve 
will be replace by the new data.  Freeze prevents this.  Up to 26 curves can be
 frozen.  
\begin{description}
\tc{Freezing Curves}
\item[(F)reeze ] This freezes the current curve 0 for the current
plotting window.  It will not be plotted in other windows. If you
change the axes from 2 to 3 dimensions and it was frozen as a 2D curve
(and {\em vice versa} ) then it will also not be plotted. It is better
to create another window to work in 3 dimensions so this is avoided.
 A parameter box
pops up that asks you for the color (linetype) as well as the key name
and the curve name.  The curve name is for easy reference and should
be a few characters.  The key name is what will be printed on the
graph if a key is present.
\item[(D)elete]  This gives you a choice of available curves to
delete.
\item[(E)dit] This lets you edit a named curve; the key, name, and
linetype can be altered.   
\item[(R)emove all]  This gets rid of all of the frozen curves in the
current window.
\item[(K)ey] This turns the key on or off.  If you turn it on, then
you can position it with the mouse on the graph.  The key consists of
a line followed by some text describing the line.  Only about 15
characters are permitted.
\item[(B)if.diag] will prompt you for a filename and then using the
current view, draw the diagram. The file must be of the same format as
is produced by the {\tt Write pts} option in the AUTO menus. (see AUTO
below.)  The diagram is colored according to whether the points are
stable/unstable fixed points or periodics. The diagram is ``frozen''
and there can only be one diagram at a time.  
\item[(C)lr. BD] clears out the current bifurcation diagram.

\end{description}
\item[a(X)es opts]  This puts up a window which allows you to tell XPP
where you want the axes to be drawn, whether you want them, and what
fontsize to make the PostScript axes labels.  
\end{description}

\item[n(U)merics] This is so important that a section is devoted to
it. See below.

\tc{File}\item[(F)ile] 
This brings up a menu with several options. Type {\tt Esc} to abort.
\begin{description}
\item[(Q)uit]  This exits XPP first asking if you are sure.
\item[(P)rt info] This prints info about the current simulation to the
console. The equations, the numerical parameters, the initial data,
boundary conditions and all parameters are printed.
\item[(W)rite set] This creates a file with all of the info about the
current numerics, etc as well as all of the currently highlighted
graphics window.  It is not meant to be readable by the user. It in
some sense saves the current state of XPP and can be read in later.
\item[(R)ead set]  This reads a set that you have previously written.
 The files are very tightly connected
 to the current ODE file so you should not load a saved file from one equation 
for a different problem.
\item[(S)ave info] This is like {\tt (P)rt info} but saves the info to
a file.  It is human readable.
\item[(C)alculator] This pops up a little window.  Type formulae in
the command 
line involving your variables and the results are displayed in the popup.  Click
 on Quit or type {\tt Esc} to exit.
\item[(A)uto] This brings up the AUTO window if you have installed
AUTO. See below for a description of this.
\tc{Edit right-hand sides} \item[(E)dit]  At long last you can edit the equations from within
XPP.  {\em Note that XPP is capable of understanding right-hand sides of up
to 256 characters.  However, the RHS editor will not accept anything
longer than about 72 characters.}  This menu item presents a
list of three options:
\begin{description}
\item[(R)HS] Edit the right-hand sides of the ODEs IDEs, and auxiliary
variables. If you are happy with the editing, then type {\tt TAB} or
click on {\tt OK.}  The program will parse the new equations and if
they are syntactically correct, alter the corresponding equation.  If
there is an error, the you will be told of the offending
right-hand-side and that will not be changed.  
\item[(F)unctions] This lets you alter any user-defined functions.  It
is otherwise the same as the above.
\item[(S)ave as] This creates an ``ODE'' file based on the current
parameter values, functions, and right-hand sides.  You will be asked
for a filename.
\end{description}
\item[C-(H)ints] spews out stuff to the terminal that represents a
skeletal C program for the right-hand sides.  Presumably, you could
use this to create faster code by replacing the file {\tt myrhs.c}
with your compiled version. Hah! I don't khow why this is even here --
Good luck!. 
\item[(B)ell on/off]  This toggles the stupid noisy bell on and off.
\item[t(I)ps] This toggles the tips on and off that appear in the
message line.
\item[(G)et par set] This loads one of the parameter sets that you
have defined in the ODE file.
\item[ (T)ranspose ]  This is not a very good place to put this but I
stuck it here just to get it into the program.  The point of this
routine is to allow one to transpose chunks of the output.  For
example, if you are solving the discretization of some spatial problem
and find a steady state, there is no way to plot the steady state as a
function of the index of the discrete system.  This routine lets you
do that.  The idea is to take something that looks like:
\begin{verbatim}
t1  x11  x21  x31 ... xm1 
t2  x12  x22  x32 ... xm2
...
tn  x1n  x2n  x3n ... xmn
\end{verbatim}
and transpose some subset of it. You are prompted for 6 items.  They
are the name of the first column you want to index, the number of
columns ({\tt ncols} and amount you want to skip across columns, {\tt
colskip} (so that {\tt colskip = 2} would be every other column.  You
must also provide the starting row {\tt j1}, the number of rows, {\tt
nrows} and the row skip, {\tt rowskip.} the The storage array is temporarily
replaced by a new array that has {\tt M=ncols} rows and {\tt nrows+1}
columns (since the data is transposed, the rows and columns are as
well; confusing ain't it).  The form of the array is:
\begin{verbatim}
1  x(i1,j1) x(i1,j2) x(i1,j3) ...
2  x(i2,j1) x(i2,j2) x(i2,j3) ...
...
M  x(iM,j1) x(iM,j2) x(iM,j3) ...
\end{verbatim}
where {\tt i2=i1+colskip, i3=i1+2*colskip, ...} and {\tt i1} is the
index corresponding to the name of the first column you
provide. Similarly, {\tt j2=j1+rowskip, ...}. As a brief example,
suppose that you solve a system of equations of the form:
\[
 x_j' = f(x_{j-1},x_j,x_{j+1},I_j)
\]
where $j=1,\dots,20.$ Click on transpose and choose {\tt x1} as the
first column, {\tt colskip=1, ncols=20} and say {\tt row1=350,
nrows=1,rowskip=1} then a new array will be produced. The first column
is the index from 1 to 20 and the second is {\tt xj(350)} where 350 is
the index and not the actual value of time. By plotting the second
column versus the first you get a ``spatial profile.'' 
\end{description}
\tc{Parameters}\item[(P)arameters]    Type the name
 of a parameter to change and enter its value. Repeat for more
 parameters.  Hit {\tt Enter} a few times to exit.   Type 
{\tt default} to get back the values when you started XPP. Use this if
 you don't want to mess with the mouse.
\tc{Erase}\item[(E)rase] erases the contents of the active window, redraws the
axes, 
deletes all 
text in the window, and sets the redraw flag to Manual.
\tc{Multiple windows}\item[(M)ake window]   The option allows
 you to create and destroy graphics windows. There are several
choices.
\begin{description}  
\item[(C)reate] makes a copy of the currently active window and makes itself active.
  You can change the graphs in this window without affecting the other windows. 

\item[(K)ill all]  Removes all but the main graphics window.

\item[(D)estroy]  This destroys the currently active window.  The main window
 cannot be destroyed.  

\item[(A)uto] turns on a flag so that the window will automatically be
redrawn when needed.

\item[(M)anual] turns off the flag and the user must restore the picture manually.

\item[(B)ottom] puts the active window on the bottom.
\end{description}
After a window is created, you can use the mouse to find the
coordinates by pressing and moving in the window.  The coordinates are
given near the top of the window.
\tc{Adding text etc}\item[(T)ext, etc] allows you to write text to the
display in a variety of sizes and in two different fonts.  You can
also add other symbols to your graph. 
\begin{description}
\item[Text] This prompts you for the text you want to add.  Then you
are asked for the size; there are five choices (0-5): 0-8pt, 1-12pt, 2-14pt,
3-18pt, 4-24pt. Next you are asked for the font: 0-TimesRoman,
1-Symbol. Note that not all X-servers will have these fonts, but the
postscript file will still draw them. Finally, place the text with the
mouse.
\item[Arrow] This lets you draw an arrow-head to indicate a direction
on a trajectory.  You will be prompted for the size, which should be
some positive number, usually less than 1.  Then you must move the the
mouse and select a direction and starting point.  Click on the
starting point and holding the mouse button down, drag the mouse to
indicate the direction of the arrow-head. Then release the
mouse-button and the arrow will be drawn.  
\item[Pointer] This is like an arrow, but draws the stem as well as
the arrow head.  It can be used to point to important features of your
graph. The prompts are like those for {\tt Arrow.}
\item[Marker] This lets you draw little markers, such as triangles,
squares, etc on the picture. When prompted to position the marker with
the mouse, you can over-ride the mouse and manually type in
coordinates if you hit the (Tab) key.
\item[markerS] This is similar to the Marker command, but allows you to
automatically mark a number of points along a computed trajectory. You
use the data browser to move the desired starting point of the list to
the top line of the browser. Then click on the (Text) (markerS)
command and choose a size and color. Then tell XPP how many markers
and how many browser lines to jump between markers. (Thus, 10 would
put a marker at every 10th data point). 

\item[Edit] This lets you edit the text, arrows, and pointers in one of
three ways:
\begin{description}
\item[Move] lets you move the object to another location without
changing any of its properties.  Choose the object with the mouse by
clicking near it.
You will then be prompted as to whether you want to move the item that
XPP selected.  If you answer {\tt yes} use the mouse to reposition it.
\item[Change] lets you change the properties: for text, the text
itself, size, and font can be change; for arrows and pointers, only
the size of the arrow head can be changed. As above, select the object
with the mouse and then edit the properties.
\item[Delete] deletes the object that you select with the mouse.
\end{description}
\item[Delete All] Deletes all the objects in the current window.
\end{description}


\tc{Equilibria}\item[(S)ing pts]  This allows you to calculate
equilibria for a discrete or continuous system.  The program also
attempts to determine stability for delay-differential equations (see
below in the numerics section.) There
 are three options.
\begin{description} 
	\item[(G)o] begins the calculation   using the values in the initial 
data box as a first guess.  Newton's method is applied.  If a value is found XPP
 tries to find the eigenvalues and asks you if you want them printed out.  If so,
 they are written to the console.  Then if there is a single real positive or
 real negative eigenvalue, the program asks you if you want the unstable or
 stable manifolds to be plotted.  Answer yes if so and they will be
 approximated. The calculation will continue until either a variable
goes out
 of bounds or you press {\tt Esc}.  If {\tt Esc} is pressed, the other
branch is
 computed. 
 The program continues to find any other invariant sets until it has gotten 
them all.  These are not stored. The program will continue integrating
beyond the current numerical parameters until you press {\tt Esc}.
Once an 
equilibrium is computed a window 
appears with info on the value of the point and its stability.  The top of 
the window tells you the number of complex eigenvalues with positive,{\tt (c+)}, 
negative {\tt (c-)}, zero {\tt (im)}  real parts and the number of
real positive
{\tt  (r+)}  and 
real negative { \tt (r-)}  eigenvalues. If the equation is a difference equation, 
then the symbols correspond the numbers of real or complex eigenvalues 
outside {\tt (+)}  the unit circle or inside { \tt (-)}
This window remains and can be iconified. 
\item[(M)ouse]  This is as above but you can specify the initial guess by 
clicking the mouse.  Only the two variables in the two-D window will reflect 
the mouse values.  This is most useful for 2D systems.
\item[(R)ange]  This allows you to find a set of equilibria over a range
 of parameters.  A parameter box will prompt you for the parameter, starting
 and ending values, number of steps.  Additionally, two other items are 
requested.  Column for stability will record the stability of the equilibrium 
point in the specified column (Use column number greater than 1).  
If you elect to shoot at each, the invariant 
manifolds will be drawn for each equilibrium computed.  The stability can be 
read as a decimal number of the form {\tt  u.s } where {\tt s}  is the
number of stable and {\tt u} 
 the number of unstable eigenvalues. So {\tt 2.03} means 3 eigenvalues
with negative real parts (or in the unit circle) and 2 with positive
real parts (outside the unit circle.)  The result of a range calculation is saved 
in the data array and replaces what ever was there.  The value of the parameter
 is in the time column, the equilibria in the remaining columns and the 
stability info in whatever column you have specified. 
 {\tt Esc} aborts one step and 
{\tt /}  aborts the whole procedure. As with the initial data/range
option, you can also make a movie.  This is useful mainly for systems
where invariant sets are to be computed.
\end{description} 
\tc{Axes}\item[(V)iew axes]  This selects one of  different types of graphs: 
a 2D box or a 3D Box; brings up a threed window; or lets you create
animations of your simulation. 
If you select the 2D curve, you will be asked for limits 
as in the window command as well as the variables to place on the axes
and the labels for the axes.   
3D is more complicated.  You will be asked for the 3 variables for the 3 axes, 
their max and min values and 4 more numbers, {\tt XLO},  etc.  XPP first scales the 
data to fit into a cube with corners (-1,-1,-1) and (1,1,1).  Rotation of this
 cube is performed and then projected into the two-D window.  {\tt XLO}, etc define
 the scales of this projection and are thus unrelated to the values of your 
data.  You are also asked for labels of the axes.
Use the {\tt (F)it } option if you don't know whats going on.
\tc{Space-time plots} (A)rray plots introduce a new window that lets the user plot many
variables at once as a function of time with color coded values.  The
point is to let one plot, e.g., an array of voltages, {\tt V1, V2,
..., VN} across the horizontal as time varies in the vertical
dimension. For example, suppose you have discretized some PDE and want
to see the evolution in space and time of the variables.  Then use
this plotting option. You will be prompted for the first column name
of the ``array'', then number of columns, the first row, then number
of rows, and the number of rows to skip. For example, if the
discretized PDE variables are {\tt u0,u1, ... ,u50}, then type in {\tt
u0} for the first element and {\tt 51} for the number of columns.
If you want rows 200
through 800 only every 4th time unit, you would put 200 for the first
row, 201 as the number of rows, and 4 as the skip value.  The {\tt Print}
button asks you for a file name and  top and bottom labels and a
render style. The render styles are:
\begin{description}
\item{-1} Grey scale
\item{0} Blue-red
\item{1} Red-Yellow-Green-Blue-Violet
\item{2} Like 1 but periodic
\end{description}
The {\tt Style} button does nothing yet. The {\tt Edit} button lets
you change ranges and arrays to plot. The {\tt Redraw} button is obvious.

Since the animation option requires learning
lots of new stuff, see section \ref{toon} for a description of the
animation language and what you can do with it.

\tc{Time plots}\item[(X)i vs t]  This chooses a certain 2D view and prompts you for the variable name. 
 The window is automatically fitted and the data plotted.  It is a
shortcut to choosing a view and windowing it.
\tc{Redrawing}\item[(R)estore] redraws the most recent data in the browser in accordance with the
 graphics parameters of the active window.
\tc{3D Parameters}\item[(3)d params]  This lets you choose rotations of the axes and
perspective  planes. 
 Play with this to see. You must be have a 3D view in the active graph
to use this. 

\tcc{3D movies!} There is another selection: {\tt Movie}.  If you choose {\tt
Yes} for this, then after you click {\tt Ok}, you will be prompted for
some additional parameters.  There are two angles you can vary, {\tt
theta} and {\tt phi}. Choose one, give an initial value, an increment,
and the number of rotations you want to perform.  XPP will then use
the {\tt Kinescope} to take successive snapshots of the screen after
performing each rotation.  You can then play these back from the
Kinescope. 
  
\tc{Boundary values}\item[(B)ndry val]  This solves boundary value problems by numerical shooting. There
 are 4 choices.
\begin{description}
\item[(S)how]  This shows the successive results of the shooting and 
erases the screen at the end and redraws the last solution.  The program
 uses the currently selected numerical integration method, the current starting
 point, {\tt T0} as the left end time and  {\tt T0+TEND} as the right
end. 
 Thus, if the
 interval of interest is {\tt (2.5,6)} then set {\tt T0=2.5} and {\tt
TEND=3.5} in the numerics menu.  
\item[(N)o show] This is as above but will not show successive solutions.
\item[(R)ange]  This allows you to range over a parameter keeping starting or
 ending values of each of the variables.  A window will appear asking you for
 the parameter, the start, end, and steps.  You will also be asked if you want 
to cycle color which means that the results of each successful solution to the
 BVP will appear in different colors. Finally the box labeled{\tt  side} tells the
 program whether to save the initial {\tt (0)} or final{\tt (1)}
values of 
the solution. 
 As the program progresses, you will see the current parameter in the info
 window under the main screen.  You can abort the current step by
pressing
 {\tt Esc} 
and the whole process by pressing {\tt /}. As in the {\tt Initialconds
Range} option, you can also choose {\tt Movie}.  Then, as before,
after each solution is computed, a snapshot is take.  Thus, you can
playback the solutions as a function of the range parameter.

\item[(P)eriodic] Periodic boundary conditions can be solved thru the
usual methods, but one then must write an addition equation for the
frequency parameter.  This option eliminates that need so that a 2-D
autonomous system need not be suspended into a 3D one.  You will be
asked for the name of the adjustable parameter for frequency.  You
will also be asked for the section variable and section.  This is an additional
condition that must be satisfied, namely, $x(0)=x_0$ where $x$ is the
section variable and $x_0$ is the section.  Type {\tt yes} if you want
the progress shown.  
\end{description}
\end{description}


\section{Numerical parameters}
When you click on the {\tt nUmerics} command in the main menu, a new
list appears.  This is the numerics menu and allows you to set all of
the numerical parameters as well as some post-processing.  Some of
these may not yet be implemented. Press {\tt Esc} or click on the {\tt
exit} to get the main menu back.

The items on this menu are:
\begin{description}
\tc{Time control}\item[(T)otal] This is the amount of time to integrate. It is called
{\tt TEND} in the documentation.  If it is
negative 
then it will be 
made positive and no data will be stored. Thus you can integrate for very long 
periods of time without being told that the storage is full.
\item[(S)tart time] This is the initial time {\tt T0} For autonomous
systems it is usually irrelevant. 
\item[tRansient] The program will integrate silently for this amount
of time before plotting output.  It is used to get rid of transients.
\item[(D)t] This is the step size used by the fixed step integrators
and the
 output step
 for Gear,CVODE, Quality RK, and Stiff algorithms.  
It is positive or negative depending on the
direction of integration.    
 \tc{Iteration control} \item[n(C)line ctrl] will prompt you for the grid size for computing
nullclines. 
\item[s(I)ng pt ctrl]  This prompts you for errors and epsilons for eigenvalues 
and equilibria calculations as well as the maximum iterates.  If you
have global flags, then you will be asked for {\tt smin} as well.

\item[n(O)ut]  This sets the number of integration steps to perform
between output 
to
 storage.  Thus, if you output every 10 steps with a {\tt Dt} of .05, XPP will yield
 output at times that are $10*.05=.5$ timesteps apart.  The advantage of this 
is that lengthier records of data can be made without losing accuracy of the 
integrator.
\item[(B)ounds]  This ses a global bound on the integrator.  If any variable exceeds
 this value in magnitude, a message appears and the integration stops.
\tc{Numerical method}\item[(M)ethod]  allows you to choose the methods of integration.
Gear, CVODE, Quality RK, and Stiff  are
 adaptive.
 Gear,CVODE, and Stiff are the best to use for stiff problems.  If you choose
the adaptive methods, you will be
 asked for error tolerance, minimum step, and maximum step size.  The discrete
 method should be used for difference equations.  If the method is
adaptive
, the
 output {\tt NOUT} is set to 1.  Also the adaptive methods generate several possible error messages
 that you may have to respond to.  These suggest how to fix the error.
 Finally, Markov processes are ignored by GEAR.
 Backward Euler prompts you for a tolerance and the maximum number of iterates
 for each step. The Volterra method is described below; it prompts for
a tolerance, maximum iterates, and a ``memory size.''  You will also
be asked if you want the convolution kernels to be re-evaluated after
any parameter is changed in either the range integration or through a
manual change in parameters.  If this flag is 1 then the kernels will
be re-computed.  The default is to not recompute them.  
\item[d(E)lay]  This sets the upper bound for the maximum delay and
three other parameters that have to do with stability of delay
equations.
 If in your delay
 equations, the delay exceeds this, a message will appear and the integration 
will stop.  Each time this is changed all previous delay data is destroyed and
 you must begin your integration anew.  Thus, it should be the first thing you
 set when solving a delay equation.  Since the storage depends on the size of
 {\tt Dt} when you change this, then the delay storage will also be
destroyed.  Delay equations require data for $t<t_0$ so that you
should edit the {\tt Delay ICs} to achieve this.  Use the fixed step
integrators and not GEAR for this.  In addition to the maximal delay,
you will be asked for {\tt AlphaMax } and {\tt OmegaMax.}  When the
program looks at stability of delay equations, it considers a
rectangular region defined by the four points in the complex plane: $
\pm i\omega_{max}, \alpha_{max}\pm i\omega_{max}.$  (See the section
on numerical methods.) Finally, you will be asked for {\tt DelayGrid}
which tells the program how many steps to take on each side of the
contour to determine stability.  Basically, the larger all these
parameters are the more accurate the stability determination.
However, if the {\tt DelayGrid} is too small relative to the size of
{\tt AlphaMax, OmegaMax} then you will also lose accuracy. 

\tc{Color coding}\item[(C)olor code]  If you have a color system, XPP can code the 
output according to
 the magnitude of the velocity or the magnitude of the Z-axis variable.  A choice
 pops up for these two options or for turning off the color.  This overrides any
 color on any other curves in your picture.  Once you choose to color code, you
 are asked to either choose max and min values or have XPP do it for you via
 optimize.  The latter will compute the max and min and set the scales accordingly.
\tc{Poincare Section}\item[(P)oincare map] This sets up parameters for Poincare sections. A choice of
 three items will appear:  the {\tt  Max/Min} option  the
 {\tt Poincare section}
 option
 and the option to turn {\tt off}  all maps. In this version, only
sections
 orthogonal
 to the coordinates are allowed, but this will change in subsequent versions.
  A parameter box will pop up and you should type in the parameters.  They are
 the variable to check and the section and the direction. That is, a point
 will be recorded if the variable crosses the section such that it is 
either positive going to negative or vice versa according to the direction
 parameter.  If the direction is set to zero, the, the point will be recorded
 from either direction.  The flag {\tt Stop on Section} instructs XPP to halt when
 the section is crossed.  Note that automatic interpolation is done.  If the 
section variable is {\tt Time, T} and the section is say {\tt T1},
then 
each time {\tt T=0 modulo T1} the point is recorded. 
 This is useful for periodically driven 
systems.  If you have opted for the {\tt Max/Min} option, then the section is 
irrelevant and the point will be recorded when a local maximum (if the 
direction is 1) minimum (direction=-1) or both (direction=0) of the variable
 is encountered.
\tc{Ruelle embedding plots}\item[R(U)elle plot]  This allows you to retard any of the axes by an integral
 number of steps.  This is useful for chaotic orbits and delayed systems. 
 Choosing a number for any of the axes will result in the variable associated
 with that axis being delayed by the number of steps inputted.  Thus if you
 plot X vs X then of course you will get a diagonal line, but if you make the 
Y-axis delayed by say 50 and the output is every .1 timesteps, then the plot 
will be X(t-5) vs X(t). This does not appear during integration of the equations
 and is available only after a computation.  You set it up and then click 
Restore from the main menu.
\tc{Stochastic stuff}\item[stoc(H)astic]  This brings up a series of items that allow you
to compute many trajectories and find their mean and variance.  It is
most useful when used with systems that are either Markovian or have
noise added to the right-hand sides.  The items are:
\begin{description}
	\item[New seed]  Use this to reseed the random number
generator.  If you use the same seed then the results will not change
from run to run.
	\item[Compute] This will put up the same dialog box as the
``Integrate'' ``Range'' choice.  Two new data sets will be created
that will compute the mean and the variance of the point by point
values of the trajectories over the number of trial runs you choose.
If the system is completely deterministic and the parameters and
initial conditions are identical for each run, then this is
superfluous.  Otherwise, the mean and variance are computed.  You can
then access these new arrays as described below.
If you fire up the sample Markov problem, choose the ``Compute''
option, and set keep the initial data constant over say 20 runs, then
you can look at the mean trajectory and its variance for each of your
variables.  
	\item[Data] This puts the results of the most recent run into
the data browser and enables plotting of them.
	\item[Mean] This puts the results of the mean value of the
most recently computed set of trials.
	\item[Variance] This does the same for the variance.
\tcc{Histogram}	\item[Histogram]  This computes a histogram for a chosen
variable and additional conditions and replaces the ``t'' column and
the first variable column with the bin values and the number per bin
respectively. You will be prompted for the number of bins, a maximum
and minimum value and the variable on which to perform the histogram.
Finally, you will be asked for additional conditions that involve the
other stored variables (not the fixed ones though.)  For example,
suppose you have run an ODE/Markov system and you want the
distribution of a continuous variable when the Markov variable is in
state 1.  Then the additional condition would be {\tt z==1} where {\tt
z} is the Markov variable.  Multiple conditions are made by using the
{\tt \&} and {\tt | } expressions.  Note that {\tt ==} is the logical
equal and is not the same as the algebraic one.  
	\item[Old Hist] brings back the most recently computed
histogram.
	\tcc{Fourier analysis} \item[Fourier] This prompts you for a
data column and the number of modes you want.  It then computes a
Fourier transform (not an FFT, since I don't want to worry about zero
padding and other matters) for the number of modes you have chosen.
The results are in the Browser.  The first column (labeled ``T'') is
the mode.  The second, the cosine component and the third, the sine
component.   
\tcc{Curve Fitting} \item[fIt curve]  This is a routine based on
Marquardt-Levenberg algorithm for nonlinear least squares fitting.
A description of the method can be found in {\it Numerical Recipes in
C.} In this implementation, one can choose parameters and initial data
to vary in an attempt to minimize the least-squares difference between
solutions to a dynamical system and data.  The data must be in a file
in which the first column contains the independent values in
increasing order.  The remaining columns contain data which are to be
fitted to solutions to a DE.  Not all the columns need be used.  When
you choose this option, a window pops up with 10 entries describing
the fit parameters. The items are:
\begin{description} 
	\item[File] This is the name of the data file. The first
column must contain the times at which the data was taken.
	\item[NCols] This contains the total number of columns in the
data file. This
includes all columns in the file, even those that you will not use.
       \item[Fitvar] This is a list separated by commas or spaces of
variables that you want to fit to the data.  These must not be Markov
variables or Auxiliary variables.  They are restricted to the items
that you define as ``Variables'' in the ODE file.  Due to laziness,
you can only have as many variables as you can type in 25 or fewer
characters.  
	\item[To Col]  This should contain a list of column numbers in
the data file associated with each of the variables you want to
fit. Thus, for example, if you want to fit ``x'' to column 5 and ``y''
to column 2, you would type ``x y'' in the ``Fitvar'' entry and ``5
2'' in the ``To Col'' entry. The number of columns in this must equal
the number of variables to be fit.
\item [Params]  These items (there are 2 of them in case you have lots
of parameters you need to vary) contain the names of parameters and
variables.  If the name is a variable, then the initial data for that
variable will be adjusted.  If it is a parameter, then the parameter
will be adjusted.  On the initial call, the current initial data and
parameter values are used. The lists of parameters and initial data
can be separated by spaces or commas.  
\item [Tolerance]  This is just a small number that tells the
algorithm when the least square error is not changing enough to be
significant.  That is if either the difference is less than ``TOL'' or
the ratio of the difference with the least square is less than ``TOL''
then the program will halt successfully.  On should not make this too
small as such differences are insignificant and a waste of CPU time.
The default of .001 seems to work well.
\item [Npts]  This is the number of points in the data file that you
want to fit to.  
\item [Epsilon]  This is used for numerical differentiation.  1e-5 is
a good value since we really don't need precise derivatives.
\item [Max iter] This is the maximum number of iterates you should use
before giving up.  
\end{description} 

On return, the program will put the best set of parameters that it has
found.  It currently is quite verbose and prints a lot of stuff to the
console.  This is mainly info about the current values of the
parameters and the least square.  
\end{description}
	


\tc{Table lookup}\item[loo(K)up ] This allows you to change the definitions of
tabulated functions by reading in a different file.  Thus, if you have
many experimental sets of data, you can read them in one by one and
integrate the equations. You are prompted for the name of a tabulated
function.  Then you give the filename to read in.  You will continue
to be prompted and can type a few carriage returns to get out. If the
table was defined as a function instead of a file, then you will be
prompted for the number of points, the limits of the range ({\tt
Xhi,Xlo}) ad finally, the formula of for the function defining the
table.  Note that it must be a function of {\tt t}.  Note that if the
function contains parameters and these are changed, it will be
automatically recomputed.


\item[bndry(V)al] prompts you for the maximum iterates, the error tolerance, and the
 deviation for the numerical Jacobian for the shooting method for solving BVPs.
\tc{Averaging}\item[(A)veraging]  This allows you to compute the adjoint and do averaging for weakly
 coupled oscillators.  To use this option, you must successfully compute
 a periodic orbit using the periodic option of the BVP menu.
Otherwise it
 is unlikely that you will obtain anything but garbage. It does not
work well with stiff systems so good luck. The menu that pops up is:
\begin{description}
	\item[(N)ew adj]  This makes the adjoint from the computed periodic
 data.  Success or failure will be noted. Separate storage is
maintained for the adjoint.  The program automatically puts the data
from the adjoint into  the browser so it can be viewed and plotted or
saved. 
\item[(A)djoint]  This will place the adjoint data in the browser,
\item[(O)rbit]  This places the periodic orbit in the browser.
\item[(M)ake H] This will prompt you for the coupling function and the result
 will be averaged and placed in 2 columns of storage, the first is the time, 
then the H function. If there are enough columns, the odd and even
parts of the averaged function will be retained. As with the adjoint,
this list is automatically placed in the browser.  
The user will be prompted for the coupling functions and 
should type in the formulas. The coupling is between two identical
units.  Say you want to couple two oscillators via a variable called
{\tt V} Then, you {\tt V} refers to the oscillator getting the input
and {\tt V'} refers to the oscillator providing it.
  For example, suppose you want to study
the behavior of two weakly diffusively coupled Fitzhugh-Nagumo
equations:
\beqa
	dv_1/dt &=& f(v_1,w_1)+\epsilon (v_2-v_1) \\
	dw_1/dt &=& g(v_1,w_1) \\
	dv_2/dt &=& f(v_2,w_2)+\epsilon (v_1-v_2) \\
	dw_2/dt &=& g(v_2,w_2) 
\eeqa
Then for the coupling you would input
\bvb
 v'-v
\end{verbatim}
\bvb
 0
\end{verbatim}
for the required coupling.  

\item[(H) function]  This places the computed H function in the browser.
\end{description}
Anytime you integrate, etc, the data will be placed back into the storage area.
\end{description}

\section{The Data Browser}
One window that pops up is the Data Browser that allows you to look at
the numbers produced by the program as well as save them to a file and
otherwise manipulate them.  It is a very primitive spread sheet.  Once
you have computed a trajectory to an equation, you can use the Data
Browser (here after, DB) to look at the data. There are some known
``bugs.'' You must sometimes grab it with the mouse and shake it to
get the data to show up (no kidding!) Also, sometimes it will not
resize properly; in this case, you must continue the integration to
the current starting point (in other words, do nothing but press C and
then <Return>.) To activate the DB, bring the window to the top and
put the pointer inside the box.  Across the top is a menu of commands
and then there follows a list of titles for the variables, time and
the auxiliary functions. Using the arrow keys, the page keys or
clicking on the appropriate commands allows you to scroll through the
data. By resizing the DB you can get more or less data.  Clicking on
Left shifts the data window to the left and Right moves it to the
right.  The time column always stays fixed.  If you do parametric or
range calculations, the range variable is kept in the time column.
Home takes you to the top of the data and End to the last row.  The
remaining commands will be described separately. The keyboard
shortcuts to invoke them are in parentheses.
\begin{description}  
\item[(F)ind] This pops up a window and asks for a variable and a
value.  It then looks through the data until it comes to the closest
value to the specified that the variable takes.  It only moves down
the file so that you can find successive values by starting at the
top.
\item[(G)et] This makes the top row of data the initial conditions
for a new run.
\item[Re(p)lace] This pops up a window asking you for a column to
replace.  Then it prompts you for a formula.  Suppose you want to
replace {\tt AUX1} with {\tt x+y-t} where {\tt x,y} are two variables.
Then type this in when prompted and the column that held {\tt AUX1}
will be replaced by the values in these columns. Any valid XPP
function or user function can be used.  Two special symbols can also
be used when applied to single variables:
\tc{Numeric integration/differentiation}\begin{description}
\item[@VARIABLE] replaces the column with the numerical derivative of
the variable.  You cannot use this within a formula, but once the
column is replaced, it can be treated as any other column.
\item[\&VARIABLE] replaces the column with the numerical integral.
Note that successive applications of the derivative and then the
integral will result in the original plus a constant.
\end{description}
\item[(U)nreplace] This undoes the most recent replacement.
\item[Fir(s)t] This marks the top row in the DB (nothing is shown)as
the start or first row for saving or restoring.
\item[Last (e)] This marks the top row as the last or end row for
saving or restoring.
\item[(T)able] This lets you save data in a tabulated format that can
then be used by XPP as a function or inputs.  You must use the {\tt
First} and {\tt Last} keys to mark the desired data you want to save.
You are then prompted for the column name, the minimum and maximum you
want your independent variable to range and the file name.  The result
is a table file of the format shown in the section on tables.
\item[(R)estore] replots the data marked by First and Last.  The
default is the entire data set
\tc{Writing output}\item[(W)rite] prompts you for a filename and writes the marked data
to a file.  The files are of the form:
\begin{verbatim}
t0 x1(t0) ... xn(t0)
t1 x1(t0) ... xn(t0)
.
.
.
tf x1(tf) ... xn(tf)
\end{verbatim}
and reflect the current contents of the DB.  They are ASCII readable.
\item[(L)oad] Will load in as much of a similarly formatted data file
as possible for graphing.

\item[(A)ddcol]  This allows you to add an additional column to the
data browser.  You are prompted for the name you want to give the
column and for the formula.  It is thus, like an auxiliary variable
with a name.  Thereafter, it will be computed along with any other
quantities that you have defined.  It is as though you had included
another auxiliary variable in your original file.
\item[(D)elcol]  This lets you delete a column.  You can only delete
columns which you have created with the {\tt Addcol} command. 

Sometimes the names will not appear or  go away on the browser.  
If you iconify it
and then open it up again, they will appear as they should.  Another
problem that can occur is that you can delete a column that is itself
referred to by a different column.  This will result in wrong answers
in the column. Thus, do not delete columns whose contents are used by
other columns.  Also, if you delete a column, its name is still known
by the internal system but it has no real value. For these reasons,
you should delete columns with caution.  If the purpose of deleting
them is to change the formula, use the right-hand-side editor,
(File-Edit) instead.



\end{description}

 
\section{Functional equations}
In the course of some research problems, I was pursuing, I ran into
some Volterra equations of the form:
\[
	u(t)=f(t)+\int_0^t K(t,s,u(s))ds
\]
that I could not convert to ODES.  (If $K$ is a convolution with a 
sum of powers and
exponentials, then it can be converted to an ODE.  Since XPP is much
more efficient with ODEs and has been thoroughly debugged with respect
to them, you should always attempt this conversion first.)  Thus, I
have added the capability to solve equations with this type of term in
them.  This has necessitated the addition of two new commands for the
ODE file and a solver for such problems. The solver is
described below in the numerical section.  Since the equation above
requires ``memory'' all the way back to $t=0$ and one often is
interested in long time behavior, XPP truncates this to:
\[
	u(t)=f(t)+\int_{\max(t-T,0)}^t K(t,s,u(s))ds
\]
where {\tt T=DT*MaxPoints}, the latter being a parameter that you
set in the numerics menu.  The default is 4000.  Thus, one
assumes that the kernel function decays for large $t$ and so the tail
will be small.  As {\tt MaxPoints} is a parameter, you can always make
it larger at the price of taking longer to evaluate right-hand sides.
Let us consider
the following equation:
\[
 u(t) = \sin(t)+\frac{1}{2}(\cos(t)-\exp(t)-t\exp(t))+\int_0^t
(t-s)\exp(-(t-s))u(s)ds
\]
whose solution is $u(t)=\sin(t).$  The following ODE file will create
this model and also add an auxiliary variable with the solution for
purposes of comparison.
\begin{verbatim}
#  voltex1.ode
u(t)=sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{(t-t')*exp(t'-t)*u}
aux utrue=sin(t)
done
\end{verbatim}
The {\tt int\{K(u,t,t')\}} construction tells XPP that this is a
Volterra integral. If your problem can be
cast as a convolution problem, considerable speedup can be obtained
since lookup tables are created.  The present example is
in fact a convolution problem, so that instead of the full
declaration, one could instead write:
\begin{verbatim}
#  voltex2.ode
u(t)=sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{t*exp(-t)#u}
aux utrue=sin(t)
done
\end{verbatim}
which convolves the first expression (of $t$ {\bf only}) with $u.$ For
this example, using a time step of .05 and integrating to 40, there is
a 3-fold speed-up using the convolution.  For more complicated
kernels, it will be more.  


 If one wants to solve, say,
\[
 u(t) = exp(-t) + \int^t_0 (t-t')^{-mu} K(t,t',u(t'))dt'
\]
the form is:
\begin{verbatim}
u(t)= exp(-t) + int[mu]{K(t,t',u}
\end{verbatim}
Note that the ``mu'' must be a number between 0 and 1.  If ``mu'' is
greater than or equal to 1, the integral is singular at 0. 

\tc{Automatic evaluation}\begin{center}
{\bf Warning}
\end{center}
If you have parameters in your definition of the kernel and you change
them, then you will have to go into the numerics menu and recompute the
kernels by calling the Method command which automatically recomputes
the kernels. Alternatively, turn on the AutoEval flag and it will be
done automatically.   

\tc{Example}I close this section with an example of a pair of coupled
oscillators in an infinite bath which has diffusion and passive decay.
The equations are:
\beqann
u(t) &=&\int_0^t k(t-s)F(u(s),v(s))ds + \int_0^t
k_d(t-s)F(u_1(s),v_1(s))ds \\
v(t) &=&\int_0^t k(t-s)G(u(s),v(s))ds + \int_0^t
k_d(t-s)G(u_1(s),v_1(s))ds \\
u_1(t) &=&\int_0^t k(t-s)F(u_1(s),v_1(s))ds + \int_0^t
k_d(t-s)F(u(s),v(s))ds \\
v_1(t) &=&\int_0^t k(t-s)G(u_1(s),v_1(s))ds + \int_0^t
k_d(t-s)G(u(s),v(s))ds 
\eeqann
where
\beqann
k(t) &=& \exp(-t)/\sqrt(\pi t) \\
k_d(t) &=& \exp(-t)\exp(-d/t)/\sqrt(\pi t)
\eeqann
and 
\beqann
F(u,v) &=& \lambda u -v - (u+qv)(u^2+v^2) \\
G(u,v) &=& \lambda v + u - (v-qu)(u^2+v^2). 
\eeqann
Note that the kernel $k$ is weakly singular and thus $\mu=.5.$

Note that in addition there is a singularity at $t=0$ for the
diffusive kernel (division by zero); this can be 
rectified by adding a small amount to the denominator. 
The XPP file is as follows
\begin{verbatim}
# lamvolt.ode 
# the four variables:   
init u=0  v=0  u1=0  v1=0  
par lam=1.5  q=0.8  d=1  u0=1  u10=0.95  
# 1/sqrt(pi)=
number spi=0.56419  
# the integral equations; since (0,0,0,0) is a rest point, I
# add a small quickly decaying transient
u(t)=u0*exp(-5*t)+spi*(int[.5]{exp(-t)#f}+int[.5]{exp(-t-d/(t+.0001))#f1})
v(t)=spi*(int[.5]{exp(-t)#g}+int[.5]{exp(-t-d/(t+.0001))#g1})
u1(t0)=u10*exp(-5*t)+spi*(int[.5]{exp(-t)#f1}+int[.5]{exp(-t-d/(t+.0001))#f})
v1(t)=spi*(int[.5]{exp(-t)#g1}+int[.5]{exp(-t-d/(t+.0001))#g})
# the four functions f,g,f1,g1
f=lam*u-v-(u*u+v*v)*(u+q*v)
g=lam*v+u-(u*u+v*v)*(v-q*u)
f1=lam*u1-v1-(u1*u1+v1*v1)*(u1+q*v1)
g1=lam*v1+u1-(u1*u1+v1*v1)*(v1-q*u1)
done
\end{verbatim}

Try it.

\section{Auto interface}

AUTO is a program that was written several years ago by Eusebius
Doedel.  It has the ability to track bifurcation curves for
steady-state and periodic systems.  The program is very powerful
particularly for following periodic orbits.  A full FORTRAN
implementation of it is available along with documentation from
Doedel.  His Email is doedel@cs.concordia.edu.  

The version supported in XPP is a subset of AUTO but allows you to do
most of the things you would nornally want to for autonomous ODEs and with 
BVPS.  In
particular, you can track fixed points, find turning points and Hopf
bifurcation points, compute two-parameter curves of turning points and
Hopf points, compute branches of periodic solutions emanating from a
Hopf point, track period-doubling bifurcations, torus bifurcations,
and two-parameter curves of fixed period orbits.  Points can be
imported into XPP as well as complete orbits.  The bifurcation
diagrams are dynamically produced and you can move around them using
the arrow keys.  Curves can be saved and reloaded for later use.
Diagrams can be saved and imported into the main XPP window. 

Click on the {\tt File Auto} menu item to bring up AUTO.

\tc{The AUTO Window}  The AUTO window consists of several parts.  The
menus are on the left. A small square in the lower left tells you
about stability.  The two windows at the bottom are respectively
information about the computed points and the hints or tips window.
This bottom window also tells you the coordinates of the main graphics
window.

\subsection{Preparation}

Before you can use AUTO, you must prepare your system for it.  You must
start your bifurcation analysis from either a  fixed point of your
model, a periodic orbit, or a solution to a boundary value problem.
  AUTO seems to work best when you start
from a steady state, but I have had success starting at periodic
orbits.  If you want to start at a steady state, find one and
integrate so that the system is at rest.  If you want to start at a
periodic orbit, then find one and make sure that the total integration
time is the ``period'' of your orbit.  This is what the AUTO interface
uses as an approximate starting period.  There are several ways to do
this; the best is to use the boundary value solver of XPP but just
plain old integration often works fine. To solve a boundary value
problem, it is necessary to find an initial set of parameters for
which you can solve the problem within XPP.  You should arrange the
``length'' of the interval to always be 1.  That is you must scale the
problem so that the domain interval of interest is $[0,1].$ You must
then compute a solution using XPP before calling AUTO.
 

For the example file {\tt lecar.ode} the parameter
of interest is {\tt iapp} and this has been set at a negative value so
that the system has a stable rest state.  The variables have been
initialized to their rest states as well.  Once you have prepared the
problem as such, you are ready to run.

Click on the ``File'' item and choose ``Auto.''  A new window will
appear which has several regions.  At the left are the 9 commands and
the ``ABORT'' key.  Below is a small window for stability information.
A large window for the bifurcation curve is labeled with the axes and
their limits.  The bottom window is for information about the points
on the diagram.  

\subsection{Choosing parameters}

The first thing you should do is tell AUTO which parameters you might
use in the bifurcation analysis. Up to 5 are allowed.  Click on
``Parameter'' and a list of 5 parameters will appear.  Type in the
names of the parameters you want to use.  For {\tt lecar.ode} use {\tt
iapp,phi,gk,vk,gna}  The default is the first 5 or fewer parameters
in your {\tt ode} file. If you have fewer than 5 parameters, only the
available ones will appear.  
\subsection{Diagram axes}

Next, you should tell AUTO the axes and the main bifurcation
parameters.  Click on ``Axes'' and 6 choices appear:
\begin{description}
\item[(H)i] This plots the maximum of the chosen variable.
\item[(N)orm]  This plots the $L_2$ norm of the solution.
\item[h(I)-lo] This plots both the max and min of the chosen variable
(convenient for periodic orbits.)
\item[(P)eriod] Plot the period versus a parameter
\item[(T)wo par] Plot the second parameter versus the primary
parameter for two-parameter continuations.
\item[(Z)oom] Use the mous to zoom in on a region.
\item[last (1) par] Use the plot parameters from the last 1-parameter
plot.
\item[last (2) par] Use plot parameters from last 2-parameter plot.
\item[(F)requency ] Plot Frequency vesus parameter.
\item[(A)verage ] Plot the average of a variable versus the parameter. 
\end{description}

After clicking, a new window pops up with the following items:
\begin{description}
\item[Y-axis] This is the variable for the y-axis of the plot.  For
two-parameter and period plots, its contents is ignored.
\item[Main Parm] This is the principal bifurcation parameter.  It must
be one of those you specified in the parameter window.  The default is
the first parameter in the parameter list.
\item[2nd Parm] This is the other parameter for two-parameter
continuations.
\item[Xmin ... Ymax] The plotting dimensions of the diagram.
\end{description}

Once you press {\tt OK} the axes will be redrawn and labeled.
For the present model, set {\tt Xmin=-.5, Ymin=-1.5, Xmax=.5,
Ymax=1.0.}
\subsection{Numerical parameters}

Next, set the NUMERICAL parameters.  When you click on this, a new
window appears with the following items:

\begin{description}
\item[Ntst] This is the number of mesh intervals for discretization
of periodic orbits.
If you are getting apparently bad results or not converging, it helps
to increase this.  For following period doubling bifurcations, it is
automatically doubled so you should reset it later.
\item[Nmax] The maximum number of steps taken along any branch.  If
you max out, make this bigger.
\item[Npr] Give complete info every {\tt Npr} steps.
\item[Ds] This is the initial step size for the bifurcation
calculation. {\em The sign of {\tt Ds} tells AUTO the direction to
change the parameter.} Since stepsize is adaptive, {\tt Ds} is just a
``suggestion.''
\item[Dsmin] The minimum stepsize (positive).
\item[Dsmax] The maximum step size.  If this is too big, AUTO will
sometimes miss important points. 
\item[Par Min] This is the left-hand limit of the diagram for the
principle parameter. The calculation will stop if the parameter is
less than this.
\item[Par Max] This is the right-hand limit of the diagram for the
principle parameter. The calculation will stop if the parameter is
greater than this.
\item[Norm Min] The lower bound for the $L_2$ norm of the solution.  If it
is less than this the calculation will stop.
\item[Norm Max] The upper bound for the $L_2$ norm of the solution.  If it
is greater than this the calculation will stop.
\end{description}

For the present model, you should set {\tt Dsmax} to be 0.05, {\tt Par
Min} to -0.45 and {\tt Par Max} to 0.45. 

\subsection{User functions}

Suppose you want to get plots at specific values of parameters or at
fixed periods of a limit cycle.  Then you can click on ``User''
which produces a menu 0-9 asking you how many points you want to keep.
Click on 0 for none or some other number.  A new window will appear
with slots for 9 items.  You can type in anything of the form:
\begin{verbatim}
	<parameter>=<value>
\end{verbatim}
or
\begin{verbatim}
	T=<value>
\end{verbatim}
AUTO will mark and save complete information for any point that
satisfies either of these criteria.  The second is used to indicate
that you want to keep a point with a particular period, e.g., {\tt
T=25} will save the any periodic orbit with period 25.  

\subsection{Running}
At this point, you are probably ready to run.  But before doing a run,
here is a hint.  You can ``save'' the diagram at this point (see below
under ``File'').  Although it is an ``empty'' diagram, all parameters
axes, and numerics are saved.  You can then reload them later on.

Click on ``Run'' to run the bifurcation.  Depending on the situation,
a number of menus can come up.  For initial exploration, there are three
choices, starting at a new steady state, periodic, or boundary value
solution. 
If you are running the
example, click on the steady-state option 
and  a nice diagram will show up and a bunch of points will move
around in the stability circle.  These indicate stability: for fixed
points, they represent exponentials of the eigenvalues; for periodics,
the Floquet multipliers. Thus those in
the circle are stable and those out of the circle are unstable.
Bifurcations occur on the circle. The outer ones are ``clipped'' so
that they will always lie in the square, thus you can keep count of
them.  

The diagram,itself, has two different lines and two different circles.
Stable fixed points are thick lines, stable periodics are solid
circles, unstable fixed points are thin lines, and unstable periodics
are open circles. Additionally, there are crosses occasionally
dispersed with numbers associated with them.  These represent
``special'' points that AUTO wants to keep.  There are several of
them:
\begin{description}
\item[EP] Endpoint of a branch
\item[LP] Limit point or turning point of a branch
\item[TR] Torus bifurcation from a periodic
\item[PD] Period doubling bifurcation
\item[UZ] User defined function
\item[MX] Failure to converge  
\item[BP] Bifurcation or branch point
\item[  ] Output every $Npr^{th}$ point.
\end{description}

\subsection{Grabbing}

You can use these special points to continue calculations with AUTO.
The ``Grab'' item lets you peruse the diagram at a leisurely pace and
to grab special points or regular points for importing into XPP or
continuing a bifurcation calculation.  Click on ``Grab'' and stuff
appears in the info window and a cross appears on the diagram.  Use
the left and right arrow keys to cruise through the diagram. The right
key goes forward and the left backward.  At the bottom, information
about the branch, the point number, the type of point, the AUTO label,
the parameters, and the period are given. The points marked by crosses
have lables and types associated with them.  The type is one of the
above.  The label corresponds to the number on the diagram.  If point
is positive, it is an unstable solution and if it is negative it is
stable. As you traverse the diagram, stability is shown in the circle.
(NOTE: On some displays, the cross does not appear; for these displays
the zoom option also does not work.  You should recompile the program
with the {\tt XORFIX} flag defined.  See the Makefile.)

You can traverse the diagram very quickly by tapping the {\tt Tab} key
which takes you the special points only.  Type {\tt Esc} to exit with
no action or type {\tt Return} to grab the point.  If it is a regular
point (i.e., not special) then the parameters and the variables will
be set to the values for that point within XPP.  You can then
integrate the equations or look at nullclines, etc.  If you grab a
special point, then you can use this as a restart point for more AUTO
calculations, such as fixed period, two-parameter studies, and
continuations. Then, you can run AUTO again.  Bifurcation diagrams are
cumulative unless you reset them in the ``File'' menu.  That is, new
stuff is continually appended to the old.  The only limit is machine
memory.

If you grab a special point and click on ``Run'' several possibilities
arise depending on the point:
\begin{description}
\item[Regular Point] Reset the diagram and begin anew.  You will be
asked first if you want to do this.
\item[Hopf Point]  
\begin{description}
	\item[Periodic] Compute the branch of periodics emanating from
the Hopf point
	\item[Extend] Continue the branch of steady states through
this point.
	\item[New Point] Restart whole calculation using this as a
starting point
	\item[Two Param] Compute a two parameter diagram of Hopf
points. 
	\end{description}
\item[Period doubling] 
	\begin{description}
	\item[Doubling] Compute the branch of period 2 solutions.
	\item[Two-param]  Compute two-parameter curve of period
doubling points.
	\end{description}
 
\item[Limit point] Compute two parameter family of limit points (fixed
points or periodic.) 
\item[Periodic point] The point is periodic so
	\begin{description}
	\item[Extend] Extend the branch
	\item[Fixed Period] Two parameter branch of fixed period
points. 
	\end{description}
	
\item[Torus point]  Compute two-parameter family of torus
bifurcations or extend the branch or compute two-parameter
fixed period.

\end{description}

{\em Before running, after a point is grabbed, be sure to set up the
correct axes and ranges for the parameters.}

\subsection{Aborting}
Any calculation can be gracefully stopped by clicking on the ``Abort''
key. This produces a new end point from which you can continue.  Note
that if there are many branches, you may have to press ``Abort''
several times. 

{\tt Clear} just erases the screen and {\tt reDraw} redraws it.  

\subsection{Saving diagrams}
{\tt File} allows you to do several things:
\begin{description}
\item[Import orbit] If the grabbed point is a special one and is a
periodic orbit, this loads the orbit into XPP for plotting.  This is
useful for unstable orbits that cant be computed by integrating
initial data.
\item[Save diagram] Writes a file for the complete diagram which you
can use later.
\item[Load Diagram] Loads a previously saved one.
\item[Postscript]  This makes a hard copy of the bifurcation diagram
\item[Reset diagram] This clears the whole thing.
\item[Write pts] This writes a file specified by the user which has 5
columns and describes the currently visible bifurcation diagram. The
first column has the coordinates of the x-axis, the second and third
columns hold the contents of the y-axis, (e.g. max and min of the
orbit). The fourth column is one of 1-4 meaning stable fixed point,
unstable fixed point, stable periodic, unstable periodic,
respectively. The fifth column is the branch number.  The main window
of XPP can import files in this format and plot them  
\end{description}

\section{Creating Animations}
\label{toon}
Many years ago, as a teenager, I used to make animated movies using
various objects like Kraft caramels (``Caramel Knowledge'') and
vegetables (``The Call of the Wild Vegetables'').  After I got my
first computer, I wanted to develop a language to automate computer
animation. As usual, things like jobs, family, etc got in the way and
besides many far better programmers have created computer assisted
animation programs.  Thus, I abandoned this idea until I recently was
simulating a simple toy as a project with an undergraduate.  I thought
it would be really cool if there were a way to pipe the output of a
solution to the differential equation into some little cartoon of the
toy. This would certainly make the visualization of the object much
more intuitive. There were immediately many scientific reasons that
would make such visualization useful as well. Watching the gaits of an
animal or the waving of cilia or the synchronization of many
oscillators could be done much better with animation
than two and three dimensional plots of the state variables. 
 
With this in mind, I have developed a simple scripting language that
allows the user to make little cartoons and show them in a dedicated
window.  The following steps are required
\begin{itemize}
\item Run the numerical simulation for however long you need it.
\item  Using a text editor create a description of the animation using
the little scripting language described below.
\item Click on the {\tt (V)iew axes} {\tt (T)oon} menu item from the
main XPP window.
\item Click on the {\tt File} item in the animation window that pops
up and give it the name of your script file.
\item If the file is OK, click on {\tt Go} item and the animation will
begin.
\end{itemize}

\subsection{The animation window}
A resizeable window will pop up when you choose the {\tt Toon} item
form the {\tt Viewaxes} submenu. This has 6 buttons on it. They are
\begin{itemize}
\item {\tt File} Click on this to load a new animation
file. Hereafter, I will call these {\tt ani} files. The usual
extension will be {\tt filename.ani}.  
\item {\tt Go} Click on here after you have loaded an ani file and it
will run.
\item {\tt Reset} This moves the animation back to the beginning
\item {\tt Fast} speeds up the animation
\item {\tt Slow} slows it down
\item {\tt Pause} pauses the animation. {\tt Go} restarts it from
where it started. This is very sluggish and it takes a long time to
respond to it. I hope to fix this someday.
\item {\tt MPEG} Clicking on this allows you to create a series of
``ppm'' files which can be gathered together into a single MPEG file
using the program {\tt mpeg\_encode.}  The resulting MPEG files are
playable by most web browsers so that you will have a permanent copy of
your animation.   {\bf WARNING!!} PPM files are very inefficient ways
of storing images but can be done quickly. A typical animation could
temporarily use 100 to 200 Mb of disk space.  Then the mpeg\_encoder
will generally compress this 100-fold and you can then delete the
``PPM'' files.    Clicking on this button results in a dialog box. You
should fill in  {\tt 1} in the {\tt Save Mpeg} entry if you want to make a
series of ``PPM'' files in order to make a movie. The {\tt Basename}
requires that you give a base filename for the series of
files. Thus,if you choose {\tt pend} then the filenames will be {\tt
pend\_0.ppm, pend\_1.ppm, } etc.  The {\tt Frameskip} entry is there in
case you do not want every single frame recorded.  So if you want
every third frame, type in {\tt 3}.
\end{itemize}
Some computer systems cannot produce the required amount of graphics
memory for a pixmap and if this happens, you cannot currently run the
animation. You will be told this and the animation window will never
appear. 

\subsection{ DASL: Dynamical Animation Scripting Language}
In order to use the animation components of XPP, you must first
describe the animation that you want to do. This is done by creating a
script with a text editor that describes the animation. You describe
the coordinates and colors of a number of simple geometric
objects. The coordinates and colors of these objects can depend on the
values of the variables (both regular and fixed, but not auxiliary) at
a given time. The animation then runs through the output of the
numerical solution and draws the objects based on that output. I will
first list all the commands and then give some examples. 

Basically,
there are two different types of objects: (i) transient and (ii)
permanent. Transient objects have coordinates that are recomputed at
every time as they are changing with the output of the
simulation. Permanent objects are computed once and are fixed through
the duration of the simulation. The objects themselves are simple
geometric figures and text that can be put together to form the
animation.

Each line in the script file consists of a command or object followed
by a list of coordinates that must be separated by semicolons. For
some objects, there are other descriptors which can be optional.
Since color is important in visualization, there are two ways a color
can be described. Either as a formula which is computed to yield a
number between 0 and 1 or as an actual color name started with the
dollar sign symbol, \$. The {\tt ani} file consists of a lines of
commands and objects which are loaded into XPP and played back on the
animation window. Here are the commands:
\begin{description}
\item{\bf dimension} xlo;ylo;xhi;yho
\item{\bf speed} delay
\item{\bf transient}
\item{\bf permanent}
\item{\bf line} x1;y1;x2;y2;color;thickness
\item{\bf rline} x1;y1;color;thickness
\item{\bf rect} x1;y1;x2;y2;color;thickness
\item{\bf frect} x1;y1;x2;y2;color
\item{\bf circ} x1;y1;rad;color;thickness
\item{\bf fcirc} x1;y1;rad;color
\item{\bf ellip} x1;y1;rx;ry;color;thickness
\item{\bf fellip} x1;y1;rx;ry;color
\item{\bf text} x1;y1;s
\item{\bf vtext} x1;y1;s;z
\item{\bf settext} size;font;color
\item{\bf end}
\end{description}

All commands can be abbreviated to their first three letters and case
is ignored. At startup the dimension of the animation window in user
coordinates is (0,0) at the bottom left and (1,1) at the top
right. Thus the point (0.5,0.5) is the center no matter what the
actual size of the window on the screen.  {\bf Color} is described by either
a floating point number between 0 and 1 with 0 corresponding to red
and 1 to violet. When described as a floating point number, it can be
a formula that depends on the variables. In all the commands, the
color is optional {\em except} {\bf settext.}  The other way of
describing color is to use names which all start with the \$
symbol. The names are:
{\bf \$WHITE, \$RED, \$REDORANGE, \$ORANGE, \$YELLOWORANGE,
                    \$YELLOW, \$YELLOWGREEN, \$GREEN, \$BLUEGREEN,
		      \$BLUE,\$PURPLE, \$BLACK}.

The {\bf transient} and {\bf permanent} declarations tell the animator
whether the coordinates have to be evaluated at every time or if they
are fixed for all time.  The default when the file is loaded is {\bf
transient.}  Thus, these are just toggles between the two different
types of objects. 

The number following the {\bf speed} declaration must be a nonnegative
integer. It tells tha animator how many milliseconds to wait between
pictures. 

The {\bf dimension} command requires 4 numbers following it. They are
the coordinates of the lower left corner and the upper right. The
defaults are (0,0) and (1,1).  

The {\bf settext} command tells the animator what size and color to
make the next text output. The size must be an integer, {\bf \{
0,1,2,3,4 \} } with 0 the smallest and 4 the biggest. The font is
either {\bf roman} or {\bf symbol.}  The color must be a named color
and not one that is evaluated. 

The remaining ten commands all put something on the screen. 
\begin{itemize}
\item {\bf line x1;y1;x2;y2;color;thickness }  draws a line from {\bf (x1,y1)} 
to {\bf (x2,y2)} in user coordinates.  These four numbers can be any
expression that involves variables and fixed variables from your
simulation. They are evaluated at each time step (unless the line is
{\bf permanent}) and this is scaled to be drawn in the window.  The
{\bf color} is optional and can either be a named color or an
expression that is to be evaluated. The {\bf thickness} is also
optional but if you want to include this, you must include the {\bf
color} as well. {\bf thickness} is any nonnegative integer and will
result in a thicker line. 

\item {\bf rline x1;y1;color;thickness } is similar to the {\bf line}
command, but a line is drawn from the endpoints of the last line drawn
to {\bf (xold+x1,yold+y1)} which becomes then new last point. All other options
are the same. This is thus a ``relative'' line.

\item {\bf rect x1;y1;x2;y2;color;thickness } draws a rectangle with
lower corner {\bf (x1,y1)} to upper corner {\bf (x2,y2)} with optional
color and thickness. 

\item {\bf frect x1;y1;x2;y2;color } draws a filled rectangle with
lower corner {\bf (x1,y1)} to upper corner {\bf (x2,y2)} with optional
color. 

\item {\bf circ x1;y1;rad;color;thick } draws a circle with radius {\bf rad}
centered at {\bf (x1,y1)} with optional
color and thickness.

\item {\bf fcirc x1;y1;rad;color } draws a filled circle with radius {\bf rad}
centered at {\bf (x1,y1)} with optional color.

\item {\bf ellip x1;y1;rx;ry;color } draws an ellipse with radii {\bf rx,ry}
centered at {\bf (x1,y1)} with optional
color and thickness. 

\item {\bf fellip x1;y1;rx;ry;color } draws a filled ellipse with
radii {\bf rx,ry}  centered at {\bf (x1,y1)} with optional
color. 

\item {\bf text x1;y1;s} draws a string {\bf s} at position {\bf
(x1,y1)} with the current color and text properties. Only the
coordinates can depend on the current values. 

\item {\bf vtext x1;y1;s;z} draws a string {\bf s} followed by
the floating point value {\bf z} at position {\bf
(x1,y1)} with the current color and text properties. Thus, you can
print out the current time or value of a variable at any given time.
\end{itemize}

{\em REMARK.}  As with lines in ODE files, it is possible to create
arrays of commands using the cobination of the {\bf [i1..i2] }
construction. Some examples are shown below. 
\subsection{Making ``Hardcopy''} 
The MPEG button allows you to save frames from the animation which can
then be compressed together to make an MPEG movie.  To do this you
must have MPEG software called {\tt mpeg\_encode.}  The versions I use
are avaliable via anonymous ftp at 
{\tt mm-ftp.cs.berkeley.edu/pub/mpeg/encode.}  Sources and binaries
are available at this site.  

Once you have {\tt mpeg\_encode} then you can make a hardcopy of the
animation by following these simple steps.  
\begin{enumerate}
\item Run the animation to see if you like it
\item Click on the MPEG button and fill it in as described above
\item Rerun the animation
\item You will now have a huge number of files called, {\em e.g.} {\tt
frame\_0.ppm, frame\_1.ppm, ... }.
\item Create an {\tt mpgfile}. I have included an example {\tt
mpgfile} in the distribution.  I also show you one below.  You need
only copy the one in the distribution and edit the lines {\tt OUTPUT}
changing the file to the desired output name and the lines after {\tt
INPUT} to the desired input filenames.  For example if you named your
PPM files {\tt bob\_0.ppm, bob\_1.ppm, ..., bob\_500 }, then change the line
after {\tt INPUT} to  {\tt bob\_*.ppm [0-500]}.  Save the {\tt
mpgfile}.
\item Type {\tt mpeg\_encode  mpgfile} and if all goes right, you will
see lots of output ending with an MPEG file that you can view
independently of XPP.
\item Delete all the {\tt *.ppm} files to give you back your disk space.
\end{enumerate}

\subsubsection{An example mpgfile}
\begin{verbatim}
# parameter file with good default values
#
# use this as a guideline for any parameters you don't really understand
# or don't care about
#
PATTERN		I
FORCE_ENCODE_LAST_FRAME
# !!! Change this line to your output 
OUTPUT		pend.mpg

BASE_FILE_FORMAT  PPM
GOP_SIZE	30
SLICES_PER_FRAME	1

PIXEL		FULL
RANGE		10
PSEARCH_ALG	LOGARITHMIC
BSEARCH_ALG	SIMPLE
IQSCALE		8
PQSCALE		10
BQSCALE		25

#IQSCALE		4
#PQSCALE		6
#BQSCALE		10


FRAME_RATE 50

REFERENCE_FRAME	ORIGINAL

#
# you really need to understand the following
#
INPUT_CONVERT	*

INPUT_DIR	.

INPUT
# change this line to your input
pend_*.ppm [0-200]
END_INPUT
\end{verbatim}

\subsection{Examples} 
I will start out with a simple pendulum example and then a bunch of
more interesting examples.
Here is the old pendulum again:
\begin{verbatim}
# damped pendulum pend.ode
dx/dt = xp
dxp/dt = (-mu*xp-m*g*sin(x))/(m*l)
pe=m*g*(1-cos(x))
ke=.5*m*l*xp^2
aux P.E.=pe
aux K.E.=ke
aux T.E=pe+ke
x(0)=2
param m=10,mu=1,g=9.8,l=1
param scale=0.008333
@ bounds=1000
done
\end{verbatim} 
I have added an initial condition and another parameter used to scale
the magnitude of the kinetic energy for a later animation. Fire up XPP
and run this simulation. Now we will create a very simple animation:
\begin{verbatim}
# pend.ani
# simple animation file for pendulum
line .5;.5;.5+.4*sin(x);.5-.4*cos(x);$BLACK;3
fcircle .5+.4*sin(x);.5-.4*cos(x);.05;$RED
end
\end{verbatim}
Notice that comments are allowed. This file is included with the
distribution as is the ODE file so you don't have to type it in. There
are only two lines of code. The first tells the animator to draw a
line from the center of the screen at {\bf (.5,.5)} to a point {\bf
( 5+.4*sin(x), .5-.4*cos(x) )}, where {\bf x} is the variable in the
ODE file for the pendulum. The line has thickness 3 and is black. The
next line of code says to draw a filled circle centered at the same
point as the line was with radius {\bf 0.05}.  This will be colored
red.  Finally, we tell the interpreter that this is the end of the
commands.
In case you haven't already done it, click on (Initialconds)
(Go) to solve the ODE. You should see a damped oscillation. Now click
on (Viewaxes) (Toon). A new window will appear with the ``test
pattern'' on the screen. In this window click on (File) and type in
{\bf pend.ani} at the prompt. XPP will tell you that two lines were
loaded successfully. Comments are ignored. Click on (Go) in the
animation window. You will see a pendulum appear and rock back and
forth. It may be somewhat jerky depending on the server and graphics
properties of your computer and graphics. You can stop it by clicking
on (Pause), speed it up by clicking on (Fast) and slow it down by
clicking on (Slow).  The reaction to mouse clicks is terribly slow, so
to test animation, I would integrate just for a short time at first
until you are satisfied. You can edit the file and reload it with the
(File) command.  

Now we consider a much more complicated animation file:
\begin{verbatim}
# pend2.ani
# fancy animation file for pendulum
PERMANENT
settext 3;rom;$PURPLE
text .25;.9;Pendulum
line 0;.5;1;.5;$BLUE;4
SPEED 10
TRANSIENT
line .5;.5;.5+.4*sin(x);.5-.4*cos(x);$BLACK;3
fcircle .5+.4*sin(x);.5-.4*cos(x);.05;1.-scale*ke
settext 1;rom;$BLACK
vtext .05;.1;t=;t
settext 1;sym;$BLACK
vtext .05;.05;q=;x
end
\end{verbatim}

The first noncomment tells the animator that what follows will be on
every frame exactly as intially defined. The text is made fairly large
and purple in Times-Roman font. The {\bf text} command puts the text a
quarter away across the screen near the top and write ``Pendulum.''
Next a thick blue line is drawn across the middle of the screen to act
as a ``tether'' for the pendulum.  We set the delay between frame to
be 10 milliseconds with the {\bf SPEED} command.  Then all the
remaining objects are to be {\bf TRANSIENT.} The first line drawn is
the arm of the pendulum. At the end, we place a filled circle, but the
color of the circle is proportional to the kinetic energy. {\bf NOTE:}
We have used the {\bf fixed variable} version of the kinetic energy,
{\bf ke}  and
not the {\bf auxiliary variable, K.E.} since the latter is not
``known'' to the internal formula compiler but the former is.
Next we set the text color and font stuff to small black roman letters
and use the {\bf vtext} command to tell the user the current
time. Finally, we set the text to symbol and plot the value of the
angle, $\theta$ which is the same key as the letter ``q.'' 

Click on the (File) button in the animator. Load the file called {\bf
pend2.ani} and run it. 

\medskip 

The next example is of a large dynamical system that represents a set
of coupled excitable cells. The ODE file is called {\tt wave.ode} and
is included in the distribution. Here it is
\begin{verbatim}
# wave.ode
# pulse wave with diffusional coupling 
param a=.1,d=.2,eps=.05,gamma=0,i=0
vv0(0)=1
vv1(0)=1
f(v)=-v+heav(v-a)
#
vv0'=i+f(vv0)-w0+d*(vv1-vv0)
vv[1..19]'=i+f(vv[j])-w[j]+d*(vv[j-1]-2*vv[j]+vv[j+1])
vv20'=i+f(vv20)-w20+d*(vv19-vv20)
#
w[0..20]'=eps*(vv[j]-gamma*w[j])
@ meth=qualrk,dt=.25,total=150,xhi=150
done
\end{verbatim}

Now we will create an animation file that plots the values of the
voltages, {\bf v0, ..., v20} 
 as beads along the vertical axis whose horizontal height is
proportional to their voltage and whose color is proportional to the
value of the recovery variables, {\bf w0, ..., w20}. Since I have run
the simulation, I know that the recovery variables are between 0 and 1
and that the voltages are between -1 and 1.  Here is the one-line
animation file for this effect. It is called {\tt wave.ani} and is
included with the distribution:
\begin{verbatim}
# wave.ani
# animated wave 
fcircle .05+.04*[0..20];.5*(vv[j]+1);.02;1-w[j]
end
\end{verbatim}
This file uses the ``array'' capabilities of the XPP parser to expand
the single line into 21 lines from 0 to 20. I just draw a filled
circle at scaled vertical coordinate and horizontal coordinate with a
small radius and colored according to the recovery variable. The
horizontal coordinate is expanded to be {\tt .05 + .04*0}, {\tt .05 +
.04*1}, etc; the vertical is {\tt .5*(vv0+1)}, {\tt .5*(vv1+1)}, etc;
and the color is {\tt 1-w0}, {1-w1}, etc. Thus this is interpreted as
a 21 line animation file.  Try it to see what it looks like. It is a
simple matter to add a vertical scale and time ticker. 

\medskip

This next example illustrates the use of the relative line
command. Here, the model is a chain of 20 oscillators representing the
phases of spinal motoneurons which control the muscles of the lamprey,
an eel-like animal.  We will also compute a cumulative ``bend'' angle
which is dependent on the phase of each oscillation. Here is the ODE
file, called {\tt lamprey.ode}
\begin{verbatim}
# example of a chain of coupled oscillators
par grad=0,phi=.1
h(u)=sin(u+phi)
par af=1,ar=1
par bend=.1
x1'=1+grad+af*h(x2-x1)
x[2..19]'=1+grad*[j]+ar*h(x[j-1]-x[j])+af*h(x[j+1]-x[j])
x20'=1+20*grad+ar*h(x9-x10)
# here is cumulative bend
an1=bend*sin(x1)
an[2..20]=bend*sin(x[j])+an[j-1]
@ bound=1000
@ total=50
done
\end{verbatim}

Fire up XPP and run this. (Note that this should be run on a torus
phase-space, but since we are only looking at the animation, it is not
important. ) Now get the animation window and load in the file {\tt
fish.ani} which looks like:
\begin{verbatim}
# fish.ani
# lamprey swimmer -- oscchain.ode
line 0;.5;.04*cos(an1);.5+.04*sin(an1);$BLACK;3
rline .04*cos(an[2..20]);.04*sin(an[j]);$BLACK;3
END
\end{verbatim}

Click on (Go) to watch it swim! I draw a series of short
line segments relative to the previous one and at an angle that is
determined by the {\bf bend} parameter in the ODE file and on the
phase of the controlling oscillator. For fun, rerun the simulation
with the parameter {\bf grad} set to 0.1.  

\medskip

The final example revisits the Lorenz equations. These equations can
be derived by looking at a simple water wheel which consists of a
circle of leaky cups with water dripping into them. (See Strogatz for
a nice derivation of the equations.) The angle of one of the cups with
respect to the viewer is found by integrating the angular velocity
which is proportional to the {\bf x} variable in the equations. Thus,
as a final example, I present an animation with 8 cups of the Lorenz
water wheel.  I have created another ODE file that has additional info
that I will use called {\tt lorenz2.ode}.  Here it is:
\begin{verbatim}
# the famous Lorenz equation set up for animated waterwheel and
# some delayed coordinates as well
init x=-7.5  y=-3.6  z=30
par r=27  s=10  b=2.66666
par c=.2  del=.1
x'=s*(-x+y)
y'=r*x-y-x*z
z'=-b*z+x*y
# x is proportional to the angular velocity so integral is angle
theta'=c*x
th[0..7]=theta+2*pi*[j]/8
# approximate the velocity vector in the butterfly coords
z1=z-del*(-b*z+x*y)
x1=x-del*(s*(-x+y))
@ dt=.025, total=40, xplot=x,yplot=y,zplot=z,axes=3d
@ xmin=-20,xmax=20,ymin=-30,ymax=30,zmin=0,zmax=50
@ xlo=-1.5,ylo=-2,xhi=1.5,yhi=2,bound=10000
done
\end{verbatim}

Fire it up with XPP and integrate it.  Then load the animation file
called {\tt lorenz.ani} and run it. Here is the file:
\begin{verbatim}
# shows the waterwheel using the integrated ang velocity
# see Strogatz book.  Use lorenz2.ode
PERMANENT 
circ .515;.515;.46;$BLACK;2
TRANSIENT
SPEED 20
frect .5+.45*sin(th[0..7]);.5+.45*cos(th[j]);.55+.45*sin(th[j]);.55+.45*cos(th[j]);$BLACK
# plotting the butterfly and a lagged version of it !!
fcirc .5+x/40;z/50;.02;$GREEN
fcirc .5+x1/40;z1/50;.02;$RED
end
\end{verbatim}

The waterwheel can be seen when running the animation. In the center
of the screen are two colored dots that represent the $(x,z)$
coordinates (green) of the attractor and the approximate velocity of
these two variables in red.










 


\section{Creating C-files for faster simulations}
For very complicated models such as discretizations of a PDE
you may want to compile the right-hand sides and run a dedicated
program for that particular model.  This is not as hard as it
seems. You must first create a library called  {\tt libxpp.a} which is
done by typing  {\tt make xpplib } after you have successfully
compiled all of XPP.  You only need to make this library once. You can
then move it to where you usually keep libraries if you want.  Now all
you do is create a C file for the right-hand sides (I will show you
below), say, it is called {\tt lorenzrhs.c} and then type:
\begin{verbatim}
gcc lorenzrhs.c -o lorenz libxpp.a -lm -lX11 
\end{verbatim}
and if all goes well, you will have an executable called ``lorenz''.
If the ODE file was called ``lorenz.ode'' then run this as follows:
\begin{verbatim}
lorenz lorenz.ode
\end{verbatim}
as XPP needs the information contained in the ODE file to tell it the
names of variables and parameters, etc.  

Since the names of the parameters are interpreted by XPP as references
to a certain array, you must communicate their values to the C
functions for your right-hand sides.  XPP provides a little utility
for this.  In the (Files) submenu. click on (c-Hints) and a bunch of
defines will be pronted to the console.  You can use these in your
program.  All of your c-files must have the following skeleton:
\begin{verbatim}
#include <math.h>

extern double constants[]; 
main(argc,argv)
 char **argv; 
 int argc;
{
 do_main(argc,argv);
 }


/* this could lead to strange animation results but 
   these will not occur if you use no fixed variables in 
   the animation script 
*/
set_fix_rhs(t,y,neq)
     int neq;
     double t,*y;
{
  return;
}


my_rhs(t,y,ydot,neq)
 double t,*y,*ydot; 
 int neq;
{

}

extra(y,t,nod,neq)
 double t,*y; 
 int nod,neq;
{

}
\end{verbatim}
In addition, you must define the parameters for the problem.  Clicking
on the (c-Hints) will essentially write this skeleton along with the
defines for the problem. Lets take the lorenz equation as an
example. Here is the ODE file:
\begin{verbatim}
# the famous Lorenz equation set up for 3d view
init x=-7.5  y=-3.6  z=30
par r=27  s=10  b=2.66666  
x'=s*(-x+y)
y'=r*x-y-x*z
z'=-b*z+x*y
@ dt=.025, total=40, xplot=x,yplot=y,zplot=z,axes=3d
@ xmin=-20,xmax=20,ymin=-30,ymax=30,zmin=0,zmax=50
@ xlo=-1.5,ylo=-2,xhi=1.5,yhi=2
done
\end{verbatim}
Run XPP with this file as the input file and click on the (File)
(c-Hints) and the following will be written to the console:
\begin{verbatim}
#include <math.h>

 extern double constants[]; 
main(argc,argv)
 char **argv; 
 int argc;
{
 do_main(argc,argv);
 }
/* defines for lorenz.ode  */ 
#define r constants[2]
#define s constants[3]
#define b constants[4]
#define X y[0]
#define XDOT ydot[0]
#define Y y[1]
#define YDOT ydot[1]
#define Z y[2]
#define ZDOT ydot[2]
my_rhs(t,y,ydot,neq)
 double t,*y,*ydot; 
 int neq;
{
  }
extra(y,t,nod,neq)
 double t,*y; 
 int nod,neq;
{
}
set_fix_rhs(t,y,neq)
     int neq;
     double t,*y;
{
}
\end{verbatim}
You just fill in the blanks as follows to produce the required C file:
\begin{verbatim}
#include <math.h>

extern double constants[];

main(argc,argv)
     char **argv;
     int argc;
{
  do_main(argc,argv);
}

/* defines for lorenz.ode  */ 
#define r constants[2]
#define s constants[3]
#define b constants[4]
#define X y[0]
#define XDOT ydot[0]
#define Y y[1]
#define YDOT ydot[1]
#define Z y[2]
#define ZDOT ydot[2]

extra(y, t,nod,neq)
     double *y,t;
     int nod,neq;
{
 return; 
}
set_fix_rhs(t,y,neq)
     int neq;
     double t,*y;
{
  return;
}
my_rhs( t,y,ydot,neq)
 double t,*y,*ydot;
 int neq;
{
 XDOT=s*(-X+Y);
YDOT=r*X-Y-X*Z;
ZDOT=-b*Z+X*Y;

}

\end{verbatim}

Now just compile this and link it with the XPP library and run it with
``lorenz.ode'' as the input and it will be a dedicated solver of the
lorenz equations.  

Here is a final example that shows you how to include user-defined
functions and auxiliary variables.  Here is the ODE file:
\begin{verbatim}
# The Morris-Lecar model as in our chapter in Koch & Segev
#  A simple membrane oscillator.  
#
params v1=-.01,v2=0.15,v3=0.1,v4=0.145,gca=1.33,phi=.333
params vk=-.7,vl=-.5,iapp=.08,gk=2.0,gl=.5,om=1
minf(v)=.5*(1+tanh((v-v1)/v2))
ninf(v)=.5*(1+tanh((v-v3)/v4))
lamn(v)= phi*cosh((v-v3)/(2*v4))
ica=gca*minf(v)*(v-1)
v'=  (iapp+gl*(vl-v)+gk*w*(vk-v)-ica)*om
w'= (lamn(v)*(ninf(v)-w))*om
aux I_ca=ica
b v-v'
b w-w'
@ TOTAL=30,DT=.05,xlo=-.6,xhi=.5,ylo=-.25,yhi=.75
@ xplot=v,yplot=w
set vvst {xplot=t,yplot=v,xlo=0,xhi=100,ylo=-.6,yhi=.5,total=100 \
	dt=.5,meth=qualrk}
done

\end{verbatim}

and here is the C-file you need:
\begin{verbatim}
#include <math.h>

extern double constants[];

main(argc,argv)
     char **argv;
     int argc;
{
  do_main(argc,argv);
}


/* defines for lecar.ode  */ 
#define v1 constants[2]
#define v2 constants[3]
#define v3 constants[4]
#define v4 constants[5]
#define gca constants[6]
#define phi constants[7]
#define vk constants[8]
#define vl constants[9]
#define iapp constants[10]
#define gk constants[11]
#define gl constants[12]
#define om constants[13]
#define V y[0]
#define VDOT ydot[0]
#define W y[1]
#define WDOT ydot[1]
#define I_CA y[2]



double minf(v)
     double v;
{
  return .5*(1+tanh((v-v1)/v2));
}

double ninf(v)
     double v;
{
  return .5*(1+tanh((v-v3)/v4));
}

double lamn(v)
     double v;
{
  return phi*cosh((v-v3)/(2*v4));
}


extra(y, t,nod,neq)
     double *y,t;
     int nod,neq;
{
  I_CA=gca*minf(V)*(V-1.0);
}

my_rhs( t,y,ydot,neq)
 double t,*y,*ydot;
 int neq;
{
VDOT=om*(iapp+gl*(vl-V)+gk*W*(vk-V)+gca*minf(V)*(1-V));
WDOT=om*lamn(V)*(ninf(V)-W);
}

\end{verbatim}

Note that the user-defined functions appear as C functions and the
auxiliary variables are almost like regular ones and are defined in
the ``extra'' function.








\section{Some comments on the numerical methods}
Most are standard.   
\begin{enumerate} 
\item The BVP solve works 
by shooting and using Newton's method.  All Jacobi matrices are computed 
numerically.  
\item The nullclines are found by dividing the window into a grid and 
evaluating the vector field at each point.  Zero contours are found
and plotted. 
\item Equilibria are found with Newton's method and eigenvalues are found by the QR
 algorithm.  Invariant sets are found by using initial data along an eigenvector
 for the corresponding eigenvalue.  The eigenvector is computed by inverse 
iteration. 
\item The Gear algorithm is out of Gear's text 
on numerical methods.
\item The two adaptive algorithms qualrk4 and stiff
 are from Numerical Recipes.
\item CVODE is based on a C-version of LSODE.  It was written by
Scott D. Cohen and Alan C. Hindmarsh, Numerical Mathematics Group,
Center for Computational Sciences and Engineering, L-316,
Lawrence Livermore National Lab, Livermore, CA 94551. email: alanh@llnl.gov.
You can get full documentation for this powerful package 
{\tt http://netlib.bell-labs.com/netlib/ode/index.html}.
  
\item Delay equations are solved by storing previous data and 
linearly interpolating from this data.  
\item Stability of delay equations is computed by a method suggested
by Tatyana Luzyanina. The linearized stability for a delay equation
results in solving a transcendental equation $f(z)=0.$ The idea is to
use the argument principle and compute the total change in the
argument of $f(z)$ as $z$ goes around a a contour $C.$ The number of
times divided by $2\pi$ tells us the number of roots of $f$ inside the
contour. Thus, XPP simply adds values of the argument of $f(z)$ at
discrete points on a large contour defined by the user and which
encloses a big chunk of the right-half plane.  Obviously the best it
can do is give sufficient conditions for instability as there could
always be roots outside the contour. But it seems to work pretty well
with modest contours except near changes in stability.   

\item Adjoints are computed for stable periodic
 orbits by integrating the negative transpose of the numerically computed
 variational equation backwards in time using backward Euler.

\item The Volterra solver uses essentially an integrator (second order)
based on the implicit product scheme described in Peter Linz's book on
Volterra equations (SIAM,1985). For ODEs implicit schemes take
considerably more time than explicit ones, but since most of the
compute time for Volterra equations is in approximating the integral,
this time penalty is minimal.  Performance is gained primarily by
taking advantage of convolution type equations.

\item Normally distributed noise is computed by the Box-Muller
transformation of uniform noise.
\item The curve-fitting is done by using a heavily customized version
of the Marquardt-Levenberg algorithm taken from Numerical Recipes in
C. 
\item The slow fourier transform is done by actually computing the
required integrals!  It is too much to ask the user to restrict the
length of her grid to be a power of two and padding with zeros sucks.
\item The algorithms in the bifurcation package are described in the
AUTO manual available from Eusebius Doedel.  I just wrote the
interface.
\end{enumerate}
 

\section{The options file}
\label{optfiles}
You can have many options files.  They are useful for initializing
XPP. However, because you can now set options from within the ODE
file, options files are probably obsolete. They have the following format:
\bvb
9x15    BIG_FONT_NAME   <-- menu fonts
fixed   SMALL_FONT_NAME <-- IC, browser, etc fonts
0       BACKGROUND (1=white,0=black)
0	IXPLT <--  X-axis variable (0=time)
1	IYPLT <--  Y-axis variable
1	IZPLT <--  Z-axis variable
0	AXES <-- type of axis (0-2d 5-3d)
1	NJMP <-- nOutput
40	NMESH <-- Nullcline mesh
4	METHOD <-- Integration method
1	TIMEPLOT <--- set to zero if one axis is not time
8000	MAXSTOR <--- maximum rows stored
20.0	TEND <-- total integration time
.05	DT <--- time step
0.0	T0 <-- start time
0.0	TRANS <-- transient
100.	BOUND <--- bounds
.0001	HMIN <-- min step for GEAR
1.0	HMAX <-- max ``   ``   ``
.00001	TOLER <-- tolerance for GEAR
0.0	DELAY <-- maximal delay
0.0	XLO  <--- 2D window sizes
20.0	XHI
-2.0	YLO
2.0	YHI
\end{verbatim}

Within an ODE file, you can call up a different options file by
typing
\bvb
option <filename>
\end{verbatim}
where {\tt <filename>} is the name of an options file.
  



\end{document}


\begin{center} {\bf Old style ODEs } \end{center}
\tc{Old style ODE files}
For  completeness, I include a description of the old style format
for ODE files.  You should not use this format anymore although it is
still supported for backward compatibility.

\bvb
<number>
change <filename>
variable <name>=<value>, ... 
...
fixed <name>,...
...
aux <name>,...
...
markov <name> <value> <nstates>
...
wiener <name>,...
...
global sign {condition} {name1=form1;name2=form2;...}
...
parameter <name>=<value>, ...
...
number <name>=<value>
...
table <name> <filename>
...
table <name> % <npts> <xlo> <xhi> <function(t)>
...
user <name> <nargs> <expression>
...
kernel <name> <mu> <expression>
...
ode <expression>
...
integral <expression>
...
ode <expression>   (for fixed and auxiliary quantities)
...
bdry <expression>
...
r <name>
{t01} {t02} ... {t0k-1}
{t10} ...
...
{tk-1,0} ...

# Comments      

done
\end{verbatim}

{\em The first line of the file contain a positive number.} This
tells XPP how much space to allocate for the value of your variables
and auxiliary functions that you might like to maintain. This can be 
up to 100 columns. This number
is the dimension of your system plus the number of extra quantities
you record.
Each of the remaining lines contains one of the above key words at its
beginning.  Since only the first letter is used, you need only use the
first letter.  
\medskip

 
The format for loading an options file is:
\bvb
c <filename>
\end{verbatim}
which loads the options file specified in {\tt<filename>}. 

The letter 'r' is used to tell XPP that you are ready to define a
tRansition table for a Markov process whose name follows.
 (The letter 'r' has no meaning
and was used since it is the first free letter I happened upon.)








\section{C Files} 
For some types of problems that have lengthy and complicated
right-hand sides, you would like to use directly compiled C code
rather than using the built in interpreter.  I have found that one
generally does not get as much speed up as merits the extra work of
recompilation for each problem.  Nevertheless, I have included an
option for doing that in XPP.  If you want to do this for a particular
ODE file, simply append the call to XPP with the option {\tt -m} which
means to make a C file.  XPP goes on as before but when you exit the
program, there will be two files called {\tt my\_cfile.c} and {\tt
my\_hfile.h.}  These files should help you in your creation of your own
file.  Basically, you should replace the call to {\tt my\_rhs()} in
 {\tt my\_rhs.c} with the contents of {\tt my\_cfile.c.} At this time,
this is incompatible with functional equations, so dont use it if you
are solving them. Nor is it compatible with tabulated stuff.


For example,
the Morris-Lecar model
\begin{verbatim}
2
variables v,w
fixed m
params v1=-.01,v2=0.15,v3=0.1,v4=0.145,gna=1.33,phi=.333
params vk=-.7,vl=-.5,iapp=.075,gk=2.0,gl=.5,om=1
user minf 1 .5*(1+tanh((arg1-v1)/v2))
user ninf 1 .5*(1+tanh((arg1-v3)/v4))
user lamn 1 phi*cosh((arg1-v3)/(2*v4))
odev  (iapp+gl*(vl-v)+gk*w*(vk-v)+gna*m*(1-v))*om
odew (lamn(v)*(ninf(v)-w))*om
odem minf(v)
b v-v'
b w-w'
done
\end{verbatim}
produces the header file:
\begin{verbatim}

#define v y__y[0]
#define w y__y[1]
double m;
#define v1 constants[1]
#define v2 constants[2]
#define v3 constants[3]
#define v4 constants[4]
#define gna constants[5]
#define phi constants[6]
#define vk constants[7]
#define vl constants[8]
#define iapp constants[9]
#define gk constants[10]
#define gl constants[11]
#define om constants[12]
double minf();
double ninf();
double lamn();

\end{verbatim}
and the C file
\begin{verbatim}

#include "my_hfile.h"
extern double constants[];
double minf(arg1)
double arg1;
{ 
 return(.5*(1+tanh((arg1-v1)/v2))
);
}

double ninf(arg1)
double arg1;
{ 
 return(.5*(1+tanh((arg1-v3)/v4))
);
}

double lamn(arg1)
double arg1;
{ 
 return(phi*cosh((arg1-v3)/(2*v4))
);
}

my_rhs(t,y__y,ydot,neq)
 double t,*y__y,*ydot;
int neq;
 { 

 do_fix();
ydot[0]= (iapp+gl*(vl-v)+gk*w*(vk-v)+gna*m*(1-v))*om;
ydot[1]=(lamn(v)*(ninf(v)-w))*om;

}
do_fix()
{
m=minf(v);

}

\end{verbatim}
Copy the file, {\tt my\_rhs.c} to some backup.
 Delete the function called {\tt my\_rhs(...)} and substitute
the contents of {\tt my\_cfile.c} putting the top declarations of {\tt
my\_cfile.c} at the beginning of {\tt my\_rhs.c}.  Then, remake XPP.
Invoke it as before by typing {\tt xpp} followed by the file name. 

\subsection{Warnings} This tool is only a guide to the production of C
code.  There are several things to watch out for:
\begin{enumerate}
\item C is case sensitive and XPP is not so make sure you are
consistent. 
\item If your parameters are variables are the same as any of the
local variables or key words in the C file {\tt my\_rhs.c} then there
will be compilation errors.  For example, don't use {\tt i.}  When in
doubt, capitalize all your parameters and variables in the ODE file
keeping the system functions such as {\tt exp} in lower case.  
\item Some XPP functions like {\tt heav} and {\tt max} are not known
to C so you must define them.  I will probably fix this later.
\item The expression {\tt y**x} or {\tt y\^x} makes no sense in C and
should be replaced by {\tt pow(y,x)}.
\end{enumerate}












\subsection{ Improved ODE parser}
The latest version of XPP/XPPAUT (1.6 and above) uses a somewhat
different syntax for ODE files.  No longer do you have to worry about
the order in which things are described, nor do you have to worry
about some names being defined before others.  The right-hand sides,
initial data, functions, etc are now much easier to input (in the
sense that they are written almost as you would in a paper.)  Many of
the same commands used in the old style format remain unchanged, but
the handling of fixed, auxiliary, and right-hand sides is much more
intuitive.  Finally, standard line continuation is allowed with the
usual UNIX character, $\backslash$.
As this new stuff has not been tested thoroughly, let me
know if you encounter problems.  Also, the old style format still
works so when in doubt use it.  

The newstyle format for the parser has the following format:

The big difference is that the name of the variable, auxiliary
quantity, fixed quantity, and Markov variable are kept with their
right-hand sides.  Thus there is no {\tt o,i,r} stuff necessary
anymore and the {\tt fixed,variable,kernel} declarations are gone.
Auxiliary variables always have user defined names and functions are
put in as you would write them. (No longer is it necessary to use {\tt
arg1,arg2,  } etc as the names of the arguments.)  For variables that
will satisfy differential equations, maps, or integral equations, you
write the equation in the most obvious fashion.  The order of
statements is unimportant (except for fixed variables which are
evaluated in the order they are defined).  Finally the standard UNIX
line continuation character $\backslash$ is recognized so you can put statements
on multiple lines.  The exception to this is the Markov transition
matrix.  Each row of the matrix must be put on the same line.














