Manual Mat Con Taug 2019

MATCONT:
Continuation toolbox for ODEs in

Matlab
W. Govaerts, Yu. A. Kuznetsov, H.G.E. Meijer,
B. Al-Hdaibat, V. De Witte, A. Dhooge, W.
Mestrom, N. Neirynck, A.M. Riet and B. Sautois
August 2019, adapted for version MatCont 7.1.
Utrecht University
The Netherlands
Belgium
The Netherlands
Contents
1 Introduction 5
1.1 Survey of functionalities supported by MatCOnt . . . . . . . . . . . . . . . 6
1.2 Testruns and how to start with Cl MatCont . . . . . . . . . . . . . . . . . 8
1.3 Computational routines in MatCont . . . . . . . . . . . . . . . . . . . . . . 11
1.4 Availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Software requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6 Disclaimer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2 Mathematical aspects of numerical continuation and handling of singulari-

ties 13
2.1 Prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.1 Pseudo-arclength continuation . . . . . . . . . . . . . . . . . . . . . . 13
2.2.2 Moore-Penrose continuation . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Stepsize control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Singularity handling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.1 Test functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.2 Multiple test functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.3 Singularity matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4.4 User location . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3 General software aspects of MatCont 18

3.1 System definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Continuation and output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3 Curve file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4 Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4.1 The options-structure . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4.2 Derivatives of the defining system of the curve . . . . . . . . . . . . . 24
3.4.3 Singularities and test functions . . . . . . . . . . . . . . . . . . . . . . 24
3.4.4 Locators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.4.5 User functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.4.6 Defaultprocessor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.7 Special processors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.8 Workspace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.9 Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4.10 Tangent search order . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4.11 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5 Failure handling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.6 General remarks on the data flow . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.7 Directories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.8 Global Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.9 Plots in MatCont . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2
4 The odefile of a dynamical system 33
4.1 Structure and construction of an odefile . . . . . . . . . . . . . . . . . . . . . 33
4.2 Handling auxiliary functions in the construction of the odefile . . . . . . . . . 38
4.3 Access to function and Jacobian values . . . . . . . . . . . . . . . . . . . . . . 38
4.4 Time integration using the odefile . . . . . . . . . . . . . . . . . . . . . . . . . 39
5 Time integration and Poincaré maps 41

5.1 Time integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.1.1 Solver output properties . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.1.2 Jacobian matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.1.3 The Steinmetz-Larter example . . . . . . . . . . . . . . . . . . . . . . 42
5.2 Poincaré section and Poincaré map . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2.1 Poincaré maps in Cl MatCont . . . . . . . . . . . . . . . . . . . . . 44
6 Equilibrium continuation 46
6.1 Mathematical definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.2 Initialization by time integration . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.3 Bifurcations and their normal form coefficients . . . . . . . . . . . . . . . . . 46
6.3.1 Branch point locator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.4 Equilibrium initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.5 Bratu example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
7 Continuation of limit cycles 53

7.1 Mathematical definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
7.2 Discretization of a limit cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
7.3 Plotting the output of a continuation of limit cycles . . . . . . . . . . . . . . 54
7.4 Initialization by time integration . . . . . . . . . . . . . . . . . . . . . . . . . 54
7.5 Bifurcations of limit cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
7.5.1 Branch Point Locator . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.5.2 Normal form coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.6 Limitcycle initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.7 Adaptive control example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
7.8 The phase response curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
8 Continuation of codim 1 bifurcations 68

8.1 Fold Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
8.1.1 Mathematical definition . . . . . . . . . . . . . . . . . . . . . . . . . . 68
8.1.2 Bifurcations along a fold curve . . . . . . . . . . . . . . . . . . . . . . 68
8.1.3 Fold initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
8.1.4 Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
8.1.5 Example: a catalytic oscillator . . . . . . . . . . . . . . . . . . . . . . 70
8.2 Hopf Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
8.2.2 Bifurcations along a Hopf curve . . . . . . . . . . . . . . . . . . . . . . 74
8.2.3 Hopf initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
8.2.4 Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.2.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3
8.3 Period Doubling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
8.3.2 Output of a continuation of period doubling bifurcation points . . . . 79
8.3.3 Bifurcations along a flip curve . . . . . . . . . . . . . . . . . . . . . . . 79
8.3.4 Period doubling initialization . . . . . . . . . . . . . . . . . . . . . . . 80
8.3.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.4 Continuation of fold bifurcation of limit cycles . . . . . . . . . . . . . . . . . 82
8.4.2 Bifurcations along a fold of cycles curve . . . . . . . . . . . . . . . . . 83
8.4.3 Fold of cycles initialization . . . . . . . . . . . . . . . . . . . . . . . . 84
8.4.4 Example: the fast Morris-Lecar equations . . . . . . . . . . . . . . . . 85
8.5 Continuation of torus bifurcation of limit cycles . . . . . . . . . . . . . . . . . 89
8.5.2 Bifurcations along a Neimark-Sacker curve . . . . . . . . . . . . . . . 90
8.5.3 Torus bifurcation initialization . . . . . . . . . . . . . . . . . . . . . . 92
8.5.4 Example: an autonomous electronic circuit . . . . . . . . . . . . . . . 92
9 Continuation of codim 2 bifurcations 98

9.1 Branch Point Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
9.1.2 Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
9.1.3 Branch Point initialization . . . . . . . . . . . . . . . . . . . . . . . . . 98
9.1.4 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
9.2 Branch Point of Cycles Continuation . . . . . . . . . . . . . . . . . . . . . . . 101
9.2.1 Mathematical Definition . . . . . . . . . . . . . . . . . . . . . . . . . . 101
9.2.2 Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
9.2.3 Branch Point of Cycles initialization . . . . . . . . . . . . . . . . . . . 102
9.2.4 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
10 Continuation of homoclinic and heteroclinic orbits 107

10.1 Homoclinic orbits: Mathematical definition . . . . . . . . . . . . . . . . . . . 107
10.1.1 Homoclinic-to-Hyperbolic-Saddle Orbits . . . . . . . . . . . . . . . . . 107
10.1.2 Homoclinic-to-Saddle-Node Orbits . . . . . . . . . . . . . . . . . . . . 109
10.2 Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
10.3 Homoclinic initialization (HHS) . . . . . . . . . . . . . . . . . . . . . . . . . . 111
10.4 Homoclinic-to-Saddle-Node initialization (HSN) . . . . . . . . . . . . . . . . . 112
10.5 CL MatCont: the MLFast example . . . . . . . . . . . . . . . . . . . . . . . 112
10.6 Heteroclinic orbits (Het) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
A Continuer example: a curve object 115
B The Brusselator example: Continuation of a solution to a boundary value

problem in a free parameter 117
4
1 Introduction
The aim of MatCont is to provide a numerical tool for the study of parameterized dynamical
systems, in particular bifurcation studies. The first MatCont - related work was done in the
master theses [35] and [32]. The package was first announced in [10] and [11]. The present
(2019) GUI was described in [34].
The present manual on MatCont is based on version 7.1 of MatCont and the runs were
tested on Matlab[27] version 9.5 (R2018b). It is meant to be a practical guide for MatCont
users without undue discussion of the details of the used algorithms. Parts of this manual
can be used as tutorials since many explicit command line runs are discussed which are also
available in the directory Testruns of MatCont, see §1.2. There are separate tutorials for
the GUI version of MatCont.
MatCont is a GUI-matlab package that builds upon a collection of routines which
can also be used independently from the command line of matlab; we refer to this use as
the command-line package CL MatCont. The aim of MatCont and CL MatCont is to
provide a continuation and bifurcation toolbox which is compatible with the standard matlab
ODE representation of differential equations. The user can easily use his/her models without
rewriting them to a specific package. The matlab programming language makes the use and
extensions of the toolbox relatively easy.
This document is structured as follows. In Chapter 2 the underlying mathematics of
continuation is treated. Section 2.4 introduces how singularities (usually, but not necessarily,
bifurcations) are handled. The toolbox specification in Chapter 3 describes the general aspects
of numerically continuing a curve in MatCont, format of the output data etcetera.
In Chapter 4 the odefile of a dynamical system is introduced. It is the handle through
which MatCont accesses the dynamical system that is being studied. The user can also
access the dynamical system directly by calling the odefile from the Matlab command line.
Chapter 5 provides more details on time integration and, more specifically, on the computation
of Poincaré sections.
The first continuation application of the toolbox, namely the initialization and continua-
tion of equilibria, with the detection of bifurcations and the computation of their normal form
coefficients, is given in Chapter 6. Chapter 7 describes the initialization and continuation of
limit cycles with the detection of bifurcations and their normal form coefficients. A special
feature is the computation of the phase response curve and its derivative.
Chapter 8 describes the continuation of codim 1 bifurcations, i.e. limit points of equilibria,
Hopf bifurcation points of equilibria, period doubling bifurcation points of limit cycles, fold
bifurcation points of limit cycles and Neimark-Sacker bifurcation points of limit cycles. Their
bifurcations are detected and normal form coefficients are computed.
Chapter 9 describes the continuation of two codim 2 bifurcations, namely branch points
of equilibria and branch points of limit cycles.
Chapter 10 deals (rather briefly) with the initialization and continuation of homoclinic
and heteroclinic orbits.
The continuation routines can also be used outside the field of parameterized dynamical
systems. We provide two examples in the appendices. A simple example of drawing a ge-
ometric curve is given in Appendix A. The continuation of a solution to a boundary value
problem in a free parameter is described in Appendix B.
5
1.1 Survey of functionalities supported by MatCOnt
Upon the development of MatCont and CL MatCont, there were multiple objectives:
• Cover as many bifurcations in ODEs as possible (all standard bifurcations with two
control parameters are now covered).
• Allow easier data exchange between programs and with matlab’s standard ODE solvers.
• Implement robust and efficient numerical methods for all computations, using minimally
extended systems where possible.
• Represent results in a form suitable for standard control, identification, and visualiza-
tion.
• Allow for easy extendibility.

In Table 1 we give a general comparison of the available features during computations for
ODEs currently supported by the most widely used software packages auto97/2000 [12],
content 1.5 [26] and MatCont/CL MatCont.
Relationships between objects of codimension 0, 1 and 2 computed by MatCont and
CL MatCont are presented in Figures 1 and 2, while the symbols and their meaning are
summarized in Tables 2 and 3, where the standard terminology is used, see [25].
codim O
0 EP LC
1 LP H LPC NS PD
2 BP CP BT ZH HH GH CPC BPC R1 R3 R4 CH LPNS PDNS R2 NSNS LPPD GPD
Figure 1: The graph of adjacency for equilibrium and limit cycle bifurcations in MatCont
An arrow in Figure 1 from O to EP or LC means that by starting time integration from

a given point we can converge to a stable equilibrium or a stable limit cycle, respectively. In
general, an arrow from an object of type A to an object of type B means that the object of type
B can be detected (either automatically or by inspecting the output) during the computation
of a curve of objects of type A. For example, the arrows from EP to H, LP, and BP mean
that we can detect H, LP and BP during the equilibrium continuation. Moreover, for each
arrow traced in the reversed direction, i.e. from B to A, there is a possibility to start the
computation of the solution of type A starting from a given object B. For example, starting
from a BT point, one can initialize the continuation of both LP and H curves. Of course, each
object of codim 0 and 1 can be continued in one or two system parameters, respectively.
6
Table 1: Supported functionalities for ODEs in auto (a), content (c) and MatCont (m).
a c m
time-integration + +
Poincaré maps +
monitoring user functions along curves computed
by continuation + + +
continuation of equilibria + + +
detection of branch points and
codim 1 bifurcations (limit and Hopf points) of equilibria + + +
computation of normal forms
for codim 1 bifurcations of equilibria + +
continuation of codim 1 bifurcations of equilibria + + +
detection of codim 2 equilibrium bifurcations
(cusp, Bogdanov-Takens, fold-Hopf,
generalized and double Hopf) + +
computation of normal forms
for codim 2 bifurcations of equilibria +
continuation of codim 2 equilibrium bifurcations
in three parameters +
continuation of limit cycles + + +
computation of phase response curves and
their derivatives +
detection of branch points and
codim 1 bifurcations (limit points, flip and
Neimark-Sacker (torus)) of cycles + + +
continuation of codim 1 bifurcations of cycles + +
branch switching at equilibrium and cycle bifurcations + + +
continuation of branch points
of equilibria and cycles +
computation of normal forms for
codim 1 bifurcations of cycles +
detection of codim 2 bifurcations of cycles +
computation of normal forms for
codim 2 bifurcations of cycles +
continuation of orbits homoclinic to equilibria + +
continuation of orbits heteroclinic to equilibria + +
7
codim
0 LC
1 HHS HSN
2 NSS NSF NFF DR* ND* TL* SH OF* IF* NCH
Figure 2: The graph of adjacency for homoclinic bifurcations in MatCont; here * stands for
S or U.
The same interpretation applies to the arrows in Figure 2, where ‘*’ stands for either S
or U, depending on whether a stable or an unstable invariant manifold is involved.
In principle, the graphs presented in Figures 1 and 2 are connected. Indeed, it is known
that curves of codim 1 homoclinic bifurcations emanate from the BT, ZH, and HH codim 2
points. The current version of MatCont fully supports, however, only one such connection:
BT → HHS.
1.2 Testruns and how to start with Cl MatCont

Cl MatCont has a large number of testruns which are collected in the directory Testruns.
One reason for that is to check that the computational routines work well in a particular
environment. But there is also a paedagogical reason.
The easiest way to start with MatCont is by using the GUI-version and going through
the tutorials. However, the command line version Cl MatCont has more functionalities
and is more flexible for advanced use.
A reasonable way to start using Cl MatCont is to first have a quick reading of the
most important (sub)sections of this manual and then go through the testruns provided in
the directory Testruns in increasing order of complexity.
As the most important (sub)sections we recommend, in that order, §1.1, §1.5, §2, §3 with
the exception of §3.6, §3.7 and §3.8, §4 with the exception of §4.3 and §4.4.
We then recommend to go through the testruns in the order suggested below and to
compare their content and their output with the discussion in the manual. Of course, it will
often be necessary to turn back to the manual for more information.
• testodefile.m, discussed in §4.3 (basic use of the system definition file)
• testtimeinteg.m and testtimeintegJacobian, discussed in §5.1.3 (time integration)
• testPoincare.m, discussed in §5.2.1 (Poincaré map)
• testmymlPRC.m, discussed in §6.2 and §7.8 (finding a stable equilibrium by time inte-
gration, continuation of equilibria, detecting a Hopf point, starting a continuation of
limit cycles, computing the phase response curve and its derivative).
8
Type of object Label
Point P
Orbit O
Equilibrium EP
Limit cycle LC
Limit Point (fold) bifurcation LP
Hopf bifurcation H
Limit Point bifurcation of cycles LPC
Neimark-Sacker (torus) bifurcation NS
Period Doubling (flip) bifurcation PD
Branch Point BP
Cusp bifurcation CP
Bogdanov-Takens bifurcation BT
Zero-Hopf bifurcation ZH
Double Hopf bifurcation HH
Generalized Hopf (Bautin) bifurcation GH
Branch Point of Cycles BPC
Cusp bifurcation of Cycles CPC
1:1 Resonance R1
1:2 Resonance R2
1:3 Resonance R3
1:4 Resonance R4
Chenciner (generalized Neimark-Sacker) bifurcation CH
Fold–Neimark-Sacker bifurcation LPNS
Flip–Neimark-Sacker bifurcation PDNS
Fold-flip LPPD
Double Neimark-Sacker NSNS
Generalized Period Doubling GPD
Table 2: Equilibrium- and cycle-related objects and their labels within the GUI
• testadaptPRC.m, discussed in §7.8 (continuation of equilibria, detection of a Hopf point,

continuation of limit cycles, starting a curve of limit cycles by starting from a limit cycle,
computation of the phase response curve and its derivative)
• testselectcycle.m, discussed in §7.4 (starting the continuation of limit cycles from

an orbit obtained by time integration).
• testbratu.m and testbratu2.m, discussed in §6.5 (user function, equilibrium contin-

uation, neutral equilibrium (H), branch point, branch switching in a branch point of
equilibria)
• testequilcataloscill.m and testLPcataloscill.m, discussed in §8.1.5 (equilibrium

continuation, detection of LP points, continuation of an LP-curve, detection of BT and
Cusp points)
• testLPHopfcataloscill.m, discussed in §8.2.5 (equilibrium continuation, detection of

Hopf points, continuation of a Hopf curve, detection of BT points, Hopf points versus
9
Type of object Label
Limit cycle LC
Homoclinic to Hyperbolic Saddle HHS
Homoclinic to Saddle-Node HSN
Neutral saddle NSS
Neutral saddle-focus NSF
Neutral Bi-Focus NFF
Shilnikov-Hopf SH
Double Real Stable leading eigenvalue DRS
Double Real Unstable leading eigenvalue DRU
Neutrally-Divergent saddle-focus (Stable) NDS
Neutrally-Divergent saddle-focus (Unstable) NDU
Three Leading eigenvalues (Stable) TLS
Three Leading eigenvalues (Unstable) TLU
Orbit-Flip with respect to the Stable manifold OFS
Orbit-Flip with respect to the Unstable manifold OFU
Inclination-Flip with respect to the Stable manifold IFS
Inclination-Flip with respect to the Unstable manifold IFU
Non-Central Homoclinic to saddle-node NCH
Table 3: Objects related to homoclinics to equilibria and their labels within the GUI
neutral equilibria, detection of GH points)
• testadapt.m, discussed in §7.7 (continuation of equilibria, detection of Hopf points)
• testadapt1.m, discussed in §7.7 (starting a curve of limit cycles from a Hopf point,
continuation of limit cycles, detection of fold and period doubling point of cycles)
• testadapt2.m, discussed in §7.7 (starting a period doubled orbit from a PD point;

detection of secondary PD points; producing the cover picture of the manual)
• testadapt3.m, discussed in §8.3.5 (starting a curve of period doubling bifurcations from

a PD point; detection of Resonance 1:2 (R2) bifurcation points)
• testEquilMLfast.m, discussed in §8.4.4 (continuation of equilibria, detection of limit

points, Hopf points and a neutral equilibrium (H))
• testLCMLfast.m, discussed in §8.4.4 (continuation of limit cycles and detection of an

LPC (limit point of cycles) point)
• testLPCMLfast.m, discussed in §8.4.4 (continuation of an LPC curve with two free

parameters).
• testtorBPC1.m and testtorBPC2.m, discussed in §8.5.4 (detection of a Neimark-Sacker

point on a curve of periodic orbits and preceding computations; in testtorBPC2.m also
the use of a user function).
• testtorBPC3.m, discussed in §8.5.4 (continuation of Neimark-Sacker bifurcations in two

free parameters)
10
• cstr1.m, cstr2 and cstr3, discussed in §9.1.4 (continuation of equilibria, detection of
LP points, continuation of an LP curve, detection of branch points with respect to two
different parameters, continuation of a BP curve with three free parameters)
• testtorBPC4.m, discussed in §9.2.4 (continuation of equilibria, detection of a Hopf point

and a branch point)
• testtorBPC5.m, discussed in §9.2.4 (continuation of limit cycles, detection of one LPC

and two BPC (branch points of cycles))
• testtorBPC6.m, discussed in §9.2.4 (branching off in a BPC, detection of a NS and a

PD bifurcation point)
• testtorBPC7.m, discussed in §9.2.4 (continuation of a BPC with three free parameters)
• testmyml.m and homoc1.m, discussed in §10.5 (computation of a curve of equilibria

from a starting point that in the first case was obtained by time integration, detection
of a Hopf point and starting a curve of limit cycles, starting a continuation of orbits
homoclinic to saddle from a limit cycle with large period)
• testdrawcurve.m, discussed in Appendix A (drawing a curve by using a continuation

algorithm)
• testbrusselator.m, discussed in Appendix B (continuation of an equilibrium solution

to a 1-dimensional PDE)
1.3 Computational routines in MatCont

MatCont is organized around the continuation of curves of 12 different types, i.e equilibrium,
limit cycle, limit point, Hopf, limit point of cycles, period doubling of cycles, Neimark-Sacker
bifurcation of cycles, branch point, branch point of cycles, orbits homoclinic to saddle, orbits
homoclinic to saddle node and heteroclinic orbits. Each of these curve types has its own
dedicated directory in MatCont, cf. §3.7.
Each curve type is related to a parameterized dynamical system and time integration
(simulation) of that system is often either desirable to confirm the results obtained from
continuation studies or needed for the initialization of the continuation curves. Therefore,
MatCont provides access to the standard Matlab integrators and two additional high-order
integrators called ode78.m and ode87.m.
Next to continuation, the main computational contribution of MatContM is in the
initializers, which come in various forms. We mention
• Some consist of a single initialization routine which can be very simple, e.g. init H EP
(starting an equilibrium curve from a Hopf point) or more complicated, e.g. init BT Hom
(starting a curve a homoclinic orbits from a Bogdanov-Takens point). We mention in
particular the quite useful initializers init LC Hom and init LC HSN in §10.3 to start
homoclinic orbits to saddle or to saddle node from a limit cycle with a long period.
• Two of the most useful initializers rely on the use of the integrator routines to find
either equilibrium points (usually stable, see §6.2) or limit cycles (usually stable, see
§7.4).
11
• Some initializers require a specific procedure, namely these for orbits homoclinic to
saddle, orbits homoclinic to saddle node and heteroclinic orbits, cf. §10.3. These have
their own dedicated directories in MatCont, cf §3.7.
Computation of the normal form coefficients of bifurcations is another important task.

Some of the involved routines are in the directories dedicated to curves on which the bifur-
cations are naturally detected. Two directories, namely MultilinearForms and LimitCycle-
Codim2 are dedicated to a more systematic handling of this type of computation.
In §5.2 we discuss the computation of Poincaré sections, based on the integrator routines.
Finally, the computation of phase response curves and their derivatives in §7.8 is a specific
feature of MatCont.
1.4 Availability
This package is freely available for download at:
http://www.sourceforge.net/
Search for ’matcont’ and then preferably go to the latest release in the directory ’matcont’.
Unzipping the downloaded file creates a directory matcont with all necessary files (see section
3.7). Download also the readme.pdf file which contains information about the release such
as where to find tutorials and the manual.
1.5 Software requirements

The present manual on MatCont is based on version 7.1 of MatCont and the runs were
tested on Matlab 9.5 (R2018b).
In principle no special matlab packages or toolboxes are necessary but the Symbolic
Toolbox will be used if it is available. We note that in the computation of normal form
coefficients derivatives up to order 5 are used and these may contain big errors when computed
by the finite difference approximations which are used if no symbolic derivatives are available.
It is important to know that the directory LimitCycle contains 7 c-files that have to be
compiled to MEX files in a platform-dependent way. The present version is expected to work
on all Windows, Unix and Mac 32-bit and 64-bit platforms when using the default c-compiler
to compile the c-files in the LimitCycle directory into MEX-files (when matcont is called for
the first time on any platform). This is no problem under Windows 32, Linux 32 and Linux
64 where the c-compiler is present by default. For Windows 64 and Mac it may be necessary
to download the compiler separately. However, for Windows32, Windows 64 and Mac64 the
MEX files are also available in the directory ”Auxiliaries” on the root of the MatCont website.
In general, compilation can depend on the Matlab version and operating system of the
computer. In case of problems, check for details provided with the MatCont release
1.6 Disclaimer
The packages MatCont and CL MatCont are freely available for non-commercial use on
an “as is” basis. In no circumstances can the authors be held liable for any deficiency,
fault or other mishappening with regard to the use or performance of MatCont and/or
CL MatCont.
For the best understanding of dynamical systems and bifurcation theory we refer to [25].
12
2 Mathematical aspects of numerical continuation and han-
dling of singularities
Consider a smooth function F : Rn+1 → Rn . We want to compute a solution curve of the
equation F (x) = 0. Numerical continuation is a technique to compute a consecutive sequence
of points which approximate the desired branch. Most continuation algorithms implement a
predictor-corrector method. The idea behind this method is to generate a sequence of points
xi , i = 1, 2, . . . along the curve, satisfying a chosen tolerance criterion: ||F (xi )|| ≤ ǫ for some
ǫ > 0 and an additional accuracy condition ||δxi || ≤ ǫ′ where ǫ′ > 0 and δxi is the last Newton
correction.
To show how the points are generated, suppose we have found a point xi on the curve.
Also suppose we have a normalized tangent vector vi at xi , i.e. Fx (xi )vi = 0, hvi , vi i = 1.
The computation of the next point xi+1 consists of 2 steps:
• prediction of a new point
• correction of the predicted point
2.1 Prediction
Suppose h > 0 which will represent a stepsize. A commonly used predictor is tangent predic-
tion:
X 0 = xi + hvi . (1)
The choice of the stepsize is discussed in section 2.3.
2.2 Correction
We assume that X 0 is close to the curve. To find the point xi+1 on the curve we use a
Newton-like procedure. Since the standard Newton iterations can only be applied to systems
with the same number of equations as unknowns, we have to append an extra scalar condition:

F (x) = 0,
(2)
g(x) = 0.
The question is how to choose the function g(x).
2.2.1 Pseudo-arclength continuation

One option for choosing g(x) is to select a hyperplane passing through X 0 that is orthogonal
to the vector vi :
g(x) = hx − X 0 , vi i . (3)
So, the Newton iteration becomes:
X k+1 = X k − Hx−1 (X k )H(X k ) (4)

F (X) Fx (X)
H(X) = , Hx (X) = . (5)
0 viT
Then one can prove that the Newton iteration for (2) will converge to a point xi+1 on the
curve from X 0 provided that the stepsize h is sufficiently small and that the curve is regular
13
V0
X0
X1
V1
X2
vi
V2
xi+1
xi vi+1
Figure 3: Moore-Penrose continuation
(rank Fx (x) = n). Having found the new point xi+1 on the curve we need to compute the
tangent vector at that point:
Fx (xi+1 )vi+1 = 0 . (6)
Furthermore the direction along the curve must be preserved: hvi , vi+1 i = 1, so we get the
(n + 1)-dimensional appended system

Fx (xi+1 ) 0
T vi+1 = . (7)
vi 1
Upon solving this system, vi+1 must be normalized.
2.2.2 Moore-Penrose continuation

CL MatCont implements a continuation method that is slightly different from the pseudo-
arclength continuation.
Definition 1 Let A be an N × (N + 1) matrix with maximal rank. Then the Moore-Penrose

inverse of A is defined by A+ = AT (AAT )−1 .
Let A be an N × (N + 1) matrix with maximal rank. Consider the following linear system
with x, v ∈ RN +1 , b ∈ RN :
Ax = b (8)
T
v x = 0 (9)
where x is a point on the curve and v its tangent vector with respect to A, i.e. Av = 0. Since
AA+ b = b and v T A+ b = hAv, (AAT )−1 bi = 0, a solution of this system is
x = A+ b. (10)
Suppose we have a predicted point X 0 using (1). We want to find the point x on the curve
which is nearest to X 0 , i.e. we are trying to solve the optimization problem:
min{||x − X 0 || | F (x) = 0} (11)

x
14
So, the system we need to solve is:
F (x) = 0 (12)
T 0
w (x − X ) = 0 (13)
where w is the tangent vector at point x. In Newton’s method this system is solved using a
linearization about X 0 . Taylor expansion about X 0 gives:
F (x) = F (X 0 ) + Fx (X 0 )(x − X 0 ) + O(||x − X 0 ||2 ) (14)
T 0 T 0 0 2
w (x − X ) = v (x − X ) + O(||x − X || ) . (15)
So when we discard the higher order terms we can see using (8) and (10) that the solution of
this system is:
x = X 0 − Fx+ (X 0 )F (X 0 ) . (16)
However, the null vector of Fx (X 0 ) is not known, therefore we approximate it by V 0 = vi ,
the tangent vector at xi . Geometrically this means we are solving F (x) = 0 in a hyperplane
perpendicular to the previous tangent vector. This is illustrated in Figure 3. In other words,
the extra function g(x) in (2) becomes:
gk (x) = hx − X k , V k i, (17)
where Fx (X k−1 )V k = 0 for k = 1, 2, . . .. Thus, the Newton iteration we are doing is:
X k+1 = X k − Hx−1 (X k , V k )H(X k , V k ) (18)
k+1 k
V = V − Hx−1 (X k , V k )R(X k , V k ) (19)

F (X) Fx (X)
H(X, V ) = , Hx (X, V ) = (20)
0 VT

Fx (X)V
R(X, V ) = . (21)
0
One can prove that under the same conditions as for the pseudo-arclength continuation, the
Newton iterations (18) and (19) converge to a point on the curve xi+1 and the corresponding
tangent vector vi+1 , respectively. In the pseudo-arclength continuation, we had to compute
a tangent vector when a new point was found. In this case however, we already compute the
tangent vectors V k at each iterate (19), so we only need to normalize the computed tangent
vectors.
2.3 Stepsize control

Stepsize control is an important issue in these algorithms. Too small stepsizes lead to unnec-
essary work being done, while too big stepsizes can lead to losing details of the curve. An
easily implementable and proven to be reliable method is convergence-dependent control.
Consider the computation of a next point using step size hi . If the computation converged,
let n denote the number of Newton iterations needed. Then the new step size hi+1 will be
selected as follows:

 hi · hdec if not converged
hi+1 = hi · hinc if converged and n < nthr (22)

hi otherwise
where hdec < 1, hinc > 1 and nthr are constants which are experimentally determined.
15
2.4 Singularity handling
This section explains the idea of singularities which can occur on a solution branch.
2.4.1 Test functions

The idea to detect singularities is to define smooth scalar functions which have regular zeros
at the singularity points. These functions are called test functions. Suppose we have a
singularity S which is detectable by a test function φ : Rn+1 → R. Also assume we have
found two consecutive points xi and xi+1 on the curve
F (x) = 0, F : Rn+1 → Rn . (23)
The singularity S will then be detected if
φ(xi )φ(xi+1 ) < 0 . (24)
Having found two points xi and xi+1 one may want to locate the point x∗ where φ(x)
vanishes. A logical solution is to solve the following system
F (x) = 0 (25)
φ(x) = 0 (26)
using Newton iterations starting at xi . However, to use this method, one should be able to
compute the derivatives of φ(x) with respect to x, which is not always easy. To avoid this
difficulty we implemented by default a one-dimensional secant method to locate φ(x) = 0
along the curve. Notice that this involves Newton corrections at each intermediate point.
2.4.2 Multiple test functions

The above is a general way to detect and locate singularities depending on one test function.
However, it may happen that it is not possible to represent a singularity with only one test
function.
Suppose we have a singularity S which depends on nt test functions. Also assume we have
found two consecutive points xi and xi+1 and all nt test functions change sign:
∀j ∈ [1, nt ] : φj (xi )φj (xi+1 ) < 0 (27)
Also assume we have found, using a one-dimensional secant method, all zeros x∗j of the test
functions. In the ideal (exact) case all these zeros will coincide:
∀j ∈ [1, nt ] : x∗ = x∗j and φj (x∗j ) = 0 (28)
Since the continuation is not exact but numerical, we cannot assume this. However, the
locations of x∗j probably will be clustered around some center point xc . In this case we will
glue the points x∗j to x∗ = xc .
A cluster will be detected if ∀i, j ∈ [1, nt ] : ||x∗i − x∗j || ≤ ǫ for some small value ǫ. In this
case we define x∗ as the mean of all located zeroes:
nt
∗ 1 X
x = x∗j (29)
nt
j=1
16
2.4.3 Singularity matrix
Until now we have discussed singularities depending only on test functions which vanish.
Suppose we have two singularities S1 and S2 , depending respectively on test functions φ1
and φ2 . Namely, assume that φ1 vanishes at both S1 and S2 , while φ2 vanishes at only S2 .
Therefore we need a possibility to represent singularities using non-vanishing test functions.
To represent all singularities we will introduce a singularity matrix (as in [26]). This matrix
is a compact way to describe the relation between the singularities and all test functions.
Suppose we are interested in ns singularities and nt test functions which are needed to
detect and locate the singularities. Then let S be the ns × nt matrix, such that:

 0 singularity i: test function j must vanish,
Sij = 1 singularity i: test function j must not vanish, (30)

otherwise singularity i: ignore test function j.
2.4.4 User location

In some cases the default location algorithm can have problems to locate a bifurcation point.
Therefore we provide the possibility to define a specific location algorithm for a particular
bifurcation. In fact, this is done for the location of branch points of equilibrium curves and
curves of limitcycles.
17
3 General software aspects of MatCont
3.1 System definition
The user defines his dynamical system in an odefile, using the framework as in the file
standard.m in the Systems directory.
In the function fun eval, the dynamical system is to be given, where the parameters
should be listed individually. Under init, the user can define some initialization parame-
ters, as the phase variable values, the timespan, etc. All phase variables and parameters
are expected to be scalar variables, not vectors or matrices. In the further functions, it is
possible to supply the symbolic derivatives of the system to various orders to increase the
speed and/or improve accuracy of the algorithm. Note that for the state variables derivatives
up to fifth order can be provided. For the parameters first order derivatives of fun eval
and of its state Jacobian can be provided. No other derivatives can be supplied through
the system definition file, since they are never used in (CL )MatCont. We recall that
the presence of symbolic derivatives, whenever necessary, typically will be stored in the sub-
fields cds.options.SymDerivate and cds.options.SymDerivativeP of cds.options. This
is done by the initializers to the curve description files. We note that 0 ≤ SymDerivative ≤ 5
denotes the order of the highest symbolically available derivative with respect to state vari-
ables; similarly 0 ≤ SymDerivativeP ≤ 2 denotes the order of the highest symbolically avail-
able derivative in which a derivative with respect to a parameter is involved. It is always
assumed that this implies the presence of the lower order derivatives.
Finally, the odefile can contain the description of any number of user functions, i.e. func-
tions that can be monitored along computed curves and whose zeros can be detected and
located. User functions can serve many purposes; an example of a sophisticated use of user
functions is given in §8.5.4.
Details and examples on the construction of the odefile of a dynamical system are given
in Chapter 4.
3.2 Continuation and output

The syntax of the continuer is:
[x,v,s,h,f] = cont(@curve, x0, v0, options);
curve is a matlab m-file where the problem is specified (cf. section 3.3). Evaluating a func-
tion by means of a function handle replaces the earlier matlab mechanism of evaluating a
function through a string containing the function name.
x0 and v0 are respectively the initial point and the tangent vector at the initial point where
the continuation starts.
options is a structure as described in section 3.4.
The arguments v0 and options can be omitted. In this case the tangent vector at x0 is
computed internally and default options are used.
The function returns a series of matrices:
x and v are the points and their tangent vectors along the curve. Each column in x and
v corresponds to a point on the curve.
18
s is an array with structures containing information about the found singularities. Its
first and last elements always refer to the first and last points of the curve, respectively, since
for convenience these are also considered ”special”.
Each element of this structure array s has the following fields:
s.index index of the singularity point in x, so s(1).index is always equal to 1 and
s(end).index is the number of computed points.
s.label label of the singularity; by convention s(1).label is ”00” and s(end).label is
”99”.
s.msg a string that may contain any information which is useful for the user,
for example the full name of the detected special point.
s.data For each special point it contains fields with additional
information, depending on the type of point, and accumulating with increasing
codimension:
• Equilibrium: s.data.v = tangent vector at the bifurcation
– Hopf point: s.data.lyapunov = first Lyapunov coefficient

– Limit point: s.data.a = normal form coefficient
∗ Bogdanov-Takens / Zero Hopf / Double Hopf / Generalized Hopf / Cusp :
s.data.c = normal form coefficient
• Limit cycle:
s.data.multipliers = multipliers at the bifurcation
s.data.timemesh = time mesh of the orbit at the bifurcation
s.data.ntst = number of test intervals
s.data.ncol = number of collocation points
s.data.parametervalues = parameter values at the bifurcation
s.data.T = period of the orbit at the bifurcation
s.data.phi = bordering vector for locating PD bifurcations
– Period-doubling point: s.data.pdcoefficient = normal form coefficient

– Limit point of cycles: s.data.lpccoefficient = normal form coefficient
– Neimark-Sacker point: s.data.nscoefficient = normal form coefficient
• Homoclinic to hyperbolic saddle / Homoclinic to saddle-node:

s.data.timemesh = time mesh of the orbit at the bifurcation
s.data.ntst = number of test intervals
s.data.ncol = number of collocation points
s.data.parametervalues = parameter values at the bifurcation
s.data.T = period of the orbit at the bifurcation
h contains some information on the continuation process. Its columns are related to the
computed points, and have the following components:
19
Stepsize Stepsize used to calculate this point (zero for initial
point and singular points)
Half the number of correction itera- For singular points this is the number of locator
tions, rounded up to the next integer iterations
User function values The values of all active user functions
Test function values The values of all active test functions
f contains different information, depending on the continuation run. For noncycle-related
continuations, the f-vector just contains the eigenvalues, if asked for. For limit cycle contin-
uations, it begins with the mesh points of the time- discretization, followed by, if they were
asked for, the PRC- and dPRC-values in all points of the periodic orbit (cf. section 7.8).
Then, if required, follow the multipliers.
It is also possible to extend the most recently computed curve with the same options (also
same number of points) as it was first computed. The syntax to extend this curve is:
[x, v, s, h, f] = cont( x, v, s, h, f, cds);
x, v, s, h and f are the results of the previous call to the continuer and cds is the global
variable that contains the curve description of the most recently computed curve (note that
this variable has to be defined as global cds in the calling command). The function returns
the same output as before, extended with the new results.
In MatCont, all curves that have been computed using a specific system are stored
in separate .mat-files, in a directory called diagram, under a subdirectory named after the
system. For example, curves of the Connor system will be kept in .mat-files under the
subdirectory Systems/Connor/diagram/. For continuation runs, each such mat-file contains
the computed x,v,s,h,f arrays, plus the cds structure and a structure related to the curve
type. Also, it contains the variables point, ctype and num. To understand their meaning,
suppose that we are computing curves of limit cycles that we start from Hopf points. The
first such computed curve then gets the name ”H LC(1)”, point stores the string ”H” and
ctype stores the string ”LC”. Furthermore, num stores the index in s of the last selected
point of the curve (the default is 1). The second curve of the same type is called ”H LC(2)”
and so on. In fact, to save storage space, only a limited number of curves of a certain type is
stored. This number can be set by the user and the default is 2. To save a computed curve
permanently, the user must change its name.
For time integration runs, cf. §5.1 (Curve Type O), the mat-file contains ctype, option,
param, point, s, t, x. Here t is the vector of time points and x is the corresponding array
of computed points. s contains data on the first and last computed points. The meaning
of point,ctype is similar to the case of continuation curves. Finally, param is the vector of
parameters of the ODE (constant during time integration) and option is a structure that
contains optional settings for time integration.
To export the computed results of a system to a different installation of MatCont one

has to copy the corresponding m-file, the mat-file and the directory of the system.
These files also contain all information needed to export the computed results to the gen-
eral matlab environment, so MatCont is really an open system.
20
MatCont also produces graphical output. 2D and 3D graphs are plotted in matlab fig-
ure windows. Such a graph can be handled as any other graph that is produced in matlab.
It can be selected using the arrow-function of the matlab figure, and the line width, line
style and colour can be altered. Markers can be set on the curve. It can be copy-pasted into
another matlab figure. In a figure, textboxes can be inserted and axes labels can be added.
Thus the user has a plethora of possibilities to combine different MatCont output graphs
into one figure.
Finally, we note that users often want to introduce new systems that are modifications of
existing systems, but with slightly different sets of state variables and/or parameters. The
best strategy to do this in MatCont is first to edit the existing system, change its name to a
new one and click ”OK” to build an m-file with a different name and no associated directory
of computed curves. Afterwards, one can edit the newly created system, make all desired
changes and click ”OK” again. For Cl MatCont see §4
3.3 Curve file

The continuer uses a special m-file where the type of solution branch is defined. This file,
further referred to as curve.m, contains the following sections:
• curve func: contains the evaluation of the right-hand side of that type of solution
branch.
• defaultprocessor: is called and executed after each point computed during a contin-
uation experiment.
• options: sets the default setting for the options-structure for this type of solution
branch (more details are given in section 3.4).
• jacobian: contains the evaluation of the jacobian of that type of solution branch.
• hessians: contains the evaluation of the hessians of that type of solution branch.
• testf: contains the test functions for detecting bifurcations along the branch.
• userf: calls the user-defined functions if there are any.
• process: is called when a bifurcation point is detected, to handle any necessary output
messages and storage.
• singmat: defines the singularity matrix of the solution branch (cf. section 2.4.3).
• locate: here specific localisation functions can be defined for bifurcations (cf. section
2.4.4)
• init: handles any special initialisations needed in the parameters or workspace.
• done: handles any special actions needed at the end of continuing the solution branch.
• adapt: this is called after every n steps, where n is user-defined. It handles any adap-
tation of parameters, subspaces, etc.
21
3.4 Options
3.4.1 The options-structure
It is possible to specify various options for the continuation run. In the continuation we use
the options structure which is initially created with contset:
options = contset;
will initialize the structure. The continuer stores the handle to the options in the variable
cds.options. Options can then be set using
options = contset(options, optionname, optionvalue);
where optionname is an option from the following list:
InitStepsize the initial stepsize (default: 0.01)
MinStepsize the minimum stepsize to compute the next point on the curve (default: 10−5 )
MaxStepsize the maximum stepsize (default: 0.1)
MaxCorrIters maximum number of correction iterations (default: 10)
MaxNewtonIters maximum number of Newton-Raphson iterations before switching to

Newton-Chords in the corrector iterations (default: 3)
MaxTestIters maximum number of iterations to locate a zero of a testfunction (default:

10)
Increment the increment to compute first order derivatives numerically (default: 10−5 )
MoorePenrose boolean indicating the use of the Moore-Penrose continuation as the Newton-
like corrector procedure (default: 1)
FunTolerance tolerance of function values: ||F (x)|| ≤ FunTolerance is the first convergence
criterium of the Newton iteration (default: 10−6 )
VarTolerance tolerance of coordinates: ||δx|| ≤ VarTolerance is the second convergence

criterium of the Newton iteration (default: 10−6 )
TestTolerance tolerance of test functions (default: 10−5 )
Singularities boolean indicating the presence of a singularity matrix (default: 0)
MaxNumPoints maximum number of points on the curve (default: 300)
Backward boolean indicating the direction of the continuation (sign of the initial tangent
vector) v0 (default: 0)
CheckClosed number of points indicating when to start to check if the curve is closed (0 =
do not check) (default: 50)
Adapt number of points after which to call the adapt-function while computing the curve
(default: 1=adapt always)
22
IgnoreSingularity vector containing indices of singularities which are to be ignored (default:
empty)
Multipliers boolean indicating the computation of the multipliers (default: 0)
Eigenvalues boolean indicating the computation of the eigenvalues (default: 0)
TSearchOrder numerical value that indicates if unit vectors are cycled in increasing order
of index (default: 1, increasing) or decreasing (set to a value different from 1), see
§3.4.10.
Userfunctions boolean indicating the presence of user functions (default: 0)
UserfunctionsInfo is an array with structures containing information about the userfunc-

tions. This structure has the following fields:
.label label of the userfunction
.name name of this particular userfunction
.state boolean indicating whether the userfunction has to be evaluated or not
PRC variable indicating the computation of the phase response curve (default: empty)
dPRC variable indicating the computation of the derivative of the phase response curve
(default: empty)
Input vector representing the input given to the system for the computation of the phase
response curve (default: 0)
This list is stored in the file contidx.m in the directory Continuer. However, options
also contains fields which are not set by the user but frozen or filled by calls to the curvefile,
namely:
SymDerivative the highest order symbolic derivative which is present (default: 0)
SymDerivativeP the highest order symbolic derivative with respect to the free parameter(s)
which is present (default: 0)
Testfunctions boolean indicating the presence of test functions (default: 0)
WorkSpace boolean indicating to initialize and clean up user variable space (default: 0)
Locators boolean vector indicating the user has provided his own locator code to locate
zeroes of test functions. Otherwise the default locator will be used (default: empty)
ActiveParams vector containing indices of the active parameter(s) (default: empty)
ActiveUParams
ActiveSParams
ActiveSParam
The last three fields are used only in the homotopy methods for the initialzation of con-
necting orbits.
Some more details follow now.
23
3.4.2 Derivatives of the defining system of the curve
In the defining system of the object that is to be continued, the derivates can be provided
that are needed for the continuation algorithm or other computations. The continuer stores
the handle to the derivatives in the variables cds.curve jacobian,cds.curve hessians.
If cds.symjac= 1, then a call to feval(cds.curve jacobian, x) must return the (n −
1) × n Jacobian matrix evaluated at point x.
If cds.symhess= 1, then a call to feval(cds.curve hessians, x) must return a 3-
2
dimensional (n − 1 × n × n) matrix H such that H(i, j, k) = ∂∂xFj ∂x
i (x)
k
.
In the present implementation in most cases cds.symhess= 0, so the ODE-file does not
provide second order derivatives, since they are not needed in the algorithms used.
3.4.3 Singularities and test functions

To detect singularities on the curve one must set the option Singularities on. Singularities
are defined using the singularity matrix, as described in section 2.4.3. The continuer stores
the handles to the singularities, the testfunctions and the processing of the singularities re-
spectively in the variables cds.curve singmat,cds.curve testf and cds.curve process.
A call to [S,L] = feval(cds.curve singmat) gets the singularity matrix S and a vector
of 2-character strings which are abbreviations of the singularities.
A call to feval(cds.curve testf, ids, x, v) then must return the evaluation of all
testfunctions, whose indices are in the integer vector ids, at x (v is the tangent vector at x).
As a second return argument it should return an array of all testfunction id’s which could
not be evaluated. If this array is not empty the stepsize will be decreased.
When a singularity is found, a call to [failed,s] = feval(cds.curve process,i,x,v,s)
will be made to process singularity i at x. This is the point where computations can be done,
like computing normal forms, eigenvalues, etc. of the singularity. These results can then be
saved in the structure s.data which can be reused for further analysis. Note that the first
and last point of the curve are also treated as singular.
3.4.4 Locators
It may be useful to have a specific locator code for locating certain singularities (cf. section
2.4.4). To use a specific locator you must set the option Locators. This is a vector in which
the index of an element corresponds to the index of a singularity. Setting the entry to 1 means
the presence of a user-defined locator. The continuer has stored the handles to the locators
in the variable cds.curve locator and will then make a call to
[x,v]=feval(cds.curve locate,i,x1,v1,x2,v2)
to locate singularity i which was detected between x1 and x2 with their corresponding tangent
vectors v1 and v2. It must return the located point and the tangent vector at that point. If
the locator was unable to find a point it should return x = [].
3.4.5 User functions

To detect zeros of userfunctions on the curve one must set the option Userfunctions on.
The continuer has stored the handles to the userfunctions cds.curve userf. First a call
to UserInfo = contget(cds.options, ’UserfunctionsInfo’, []) is made to get infor-
mation on the userfunctions. A call to feval(cds.curve userf, UserInfo, ids, x, v)
24
then must return the evaluation of all userfunctions ids, whose information is in the structure
UserInfo, at x (v is the tangent vector at x). As a second return argument it should return
an array of all user function id’s which could not be evaluated. If this array is not empty the
stepsize will be decreased.
A special point on a bifurcation curve that is specified by a user function has a structure
as follows:
s.index index of the detected singular point defined by the user function.
s.label a string that is in UserInfo.label, label of the singularity.
s.msg a string that is set in UserInfo.name.
s.data an empty tangent vector or values of the user functions in the singular point.
When a change of sign of a userfunction is detected, the userfunction i is processed at x.
This is the point where the results (values of the userfunction) can be saved in the structure
s.data which can be reused for further analysis.
3.4.6 Defaultprocessor
In many cases it is useful to do some general computations for every calculated point on
the curve. The results of these computations can then be used by for example the test-
functions. The continuer has stored the handle to the defaultprocessor in the variable
cds.curve defaultprocessor.
The defaultprocessor is called as
[failed,f,s] = feval(cds.curve defaultprocressor,x,v,s).
x and v are the point on the curve and it’s tangent vector. The argument s is only supplied
if the point is a singular point, in that case the defaultprocessor may also add some data to
the s.data field. If for some reason the default processor fails it should set failed to 1. This
will result in a reduction of the stepsize and a retry which should solve the problem. Any
information that is to be preserved, should be put in f. f must be a column vector and must
be of equal size for every call to the default processor.
3.4.7 Special processors

After a singular point has been detected and located a singular point data structure will
be created and initialized as described in section 3.2. If there are some special data (like
eigenvalues) which may be of interest for a particular singular point then a call to [failed,s]
= feval(cds.curve process,i,x,v,s) should store this data in the s.data field. Here i
indicates which singularity was detected and x and v are the point and tangent vector where
this singularity was detected.
3.4.8 Workspace
During the computation of a curve it is sometimes necessary to introduce variables and do
additional computations that are common to all points of the curve. The continuer has stored
the handle to the initialization and cleaning of the workspace in the variables cds.curve init
and cds.curve done. These can be relegated to a call of the type
feval(cds.curve_init,x,v).
This option has to be provided only if the variable WorkSpace in cds.options is switched
on. In this case a call
25
feval(cds.curve_done,x,v)
must clear the workspace. Variables in the workspace must be set global.
3.4.9 Adaptation
It is possible to adapt the problem while generating the curve. If Adapt has a value, say 5,
then after 5 computed points a call to [reeval,x,v]=feval(cds.curve adapt,x,v) will be
made where the user can program to change the system.
For some applications it is useful to change or modify the used test functions while com-
puting the curve (like in bordering techniques). In order to preserve the correct signs of the
test functions it is sometimes necessary to reevaluate the test functions after adaptation. To
do this reeval should be one, otherwise zero. The return variables x and v should be the
updated x and v which may have changed because of the changes made to the system.
3.4.10 Tangent search order

To start a continuation, an initial point x0 and a tangent vector v0 are needed in gen-
eral. Often, only x0 is available. In this case, MatCont successively tries all unit vec-
tors as candidate tangent vectors. By default, this is done in increasing order of index
(cds.options.TSearchOrder = 1). If cds.options.TsearchOrder is set to a value different from
1 then the cycling is done in decreasing order of index. See §7.4 for an example.
In cases where the number of continuation variables is large (e.g. when computing limit
cycles) the choice of cds.options.TSearchOrder can substantially change the speed of the
computation.
3.4.11 Summary
In the following table we list calls that can be made to the continuation curve description cds
and which options are involved.
26
Syntax of call What it should do (options involved)
feval(cds.curve func,x) return F (x)
feval(cds.curve options) return option vector
feval(cds.curve jacobian,x) return Jacobian at x (SymDerivative≥ 1)
feval(cds.curve hessians,x) return Hessians at x (SymDerivative≥ 2)
feval(cds.curve init,x,v) initialize user variable space (WorkSpace)
feval(cds.curve done) destroy user variable space (WorkSpace)
feval(cds.curve defaultprocessor,x,v,s) initialize data for testfunctions and
set some general singularity data
feval(cds.curve testf,ids,x,v) return evaluation of testfunctions ids at x
(Singularities)
feval(cds.curve locate,i,x1,x2,v1,v2) return located singularity and tangent
vector(Locators)
feval(cds.curve userf,UserInfo,ids,x,v) return evaluation of userfunctions ids with
UserInfo at x (Userfunctions)
feval(cds.curve singmat) return singularity matrix (Singularities)
feval(cds.curve process,i,x) run processor code of singularity i
at x(Singularities)
feval(cds.curve adapt,x,v) run adaptation code of problem (Adapt)
3.5 Failure handling

During the continuation numerical problems may arise. For example a linear system in the
Newton corrections could be ill conditioned. In such cases the continuer checks for the last
warning issued by matlab using lastwarn() and decreases the step size along the curve.
The same mechanism is used when a test function or a user function cannot be computed.
3.6 General remarks on the data flow

At this point we have discussed two components of a continuation process, the continuer
itself and the curve definition. In Figure 4 the complete structure is visualized. The arrows
show the flow of information between the objects, where an arrow from object A to object
B indicates that information present in A is sent to B, typically by a call from B to A. The
information is sometimes passed via a function call but in many cases via a global structure.
Global structures are discussed in §3.8.
As one can see, two extra components are included: the curve initializer and some external
ODE file.
Continuation of curves with complicated curve definitions often needs to be initialized.
Since the continuer is called only with the start point x0 and an options structure (and some-
times, but not always v0 ) there must be some way to initialize other parameters. Calling
an initializer from a GUI or command prompt solves this problem. The interaction between
initializer and continuer is “invisible” since it passes through a global structure called cds
(continuation descriptor structure), see §3.8. One important field is cds.symjac which in-
forms the continuer whether or not the curve definition file includes the Jacobian of the curve
definition function. See also the note at the end of §A.
The standard matlab odeget and odeset only support Jacobian matrices coded in the
27
ODE FILE
CURVE DEFINITION
CONTINUER CURVE INITIALIZER
MATLAB PROMPT / GUI
Figure 4: Structure of continuation process
ode-file. However, we do need the derivatives with respect to the parameters. It is also useful
to have higher-order symbolic derivatives available.
To overcome this problem, the package contains new versions of odeget and odeset which
support Jacobians with respect to parameters and higher-order derivatives. The new routines
are compatible with the ones provided by matlab.
To include the Jacobian with respect to parameters, the option JacobianP should contain
the handle of the subfunction jacobianp @jacobianp. A call to feval(@jacobianp, 0, x,
p1, p2, ...) should then return the Jacobian with respect to to parameter p1 , p2 , . . . .
To include Hessians in the ode-file the option Hessians should contain the handle of the
subfunction hessians @hessians. The software then assumes that a call to feval(@hessians,
0, x, p1, p2, ...) will return all Hessians in the same way as mentioned above. Setting
the option to [] indicates that there are no Hessians available from the ode-file (default
behaviour).
To include Hessians with respect to parameters in your ode-file the option HessiansP
should contain the handle of the subfunction hessiansp @hessiansp. The software then as-
sumes that a call to feval(@hessiansp, 0, x, p1, p2, ...) will return all Hessians with
respect to parameters in the same way as mentioned above. Setting the option to [] indicates
that there are no Hessians with respect to parameters available from the ode-file (default
behaviour).
To include the third order derivatives in your ode-file the option Der3 should contain the
handle of the subfunction der3 @der3. The software then assumes that a call to feval(@der3,
0, x, p1, p2, ...) will return all third order derivatives in the same way as mentioned
above. Setting the option to [] indicates that they are not available from the ode-file (default
28
behaviour)
Der4 and Der5 are values indicating the 4th and 5th order symbolic derivative, available
in the ode-file.
3.7 Directories
The files of the toolbox are organized in the following 23 directories
• BranchPoint
Here are all files stored needed to do a branch point continuation. This includes the
initializers and a branch point curve definition file. BranchPoint and BranchPointCycle
are the only continuation curve types with three free system parameters in MatCont.
• BranchPointCycle
Here are all files stored needed to do a branch point of cycles continuation. This inluded
the initializers and branch point of cycles definition files. BranchPoint and Branch-
PointCycle are the only continuation curve types with three free system parameters in
MatCont.
• Continuer
Here are all the main files stored for the continuer, which are needed to calculate and
plot any curve.
• Equilibrium
Here are all files stored needed to do an equilibrium continuation. This includes the
initializers and the equilibrium curve definition file.
• GUI
This directory contains all GUI-related files. However, some files in this directory are
also used in Cl MatCont.
• Help
Contains the help-files.
• Heteroclinc
Here are all files stored needed to do a continuation of heteroclinic orbits. This includes
the initializers and the curve definition file.
• Homoclinic
Here are all files stored needed to do a homoclinic-to-hyperbolic-saddle continuation.
This includes the initializers and the curve definition file.
• HomoclinicSaddleNode
Here are all files stored needed to do a homoclinic-to-saddle-node continuation. This
includes the initializers and the curve definition file.
• HomtopyHet
Here are all files stored needed to find a heteroclinic connection by the homotopy
method.
29
• HomotopySaddle
Here are all files stored needed to find an orbit homoclinic to saddle by the homotopy
method.
• HomotopySaddlenode
Here are all files stored needed to find an orbit homoclinic to saddle-node by the homo-
topy method.
• Hopf
Here are all files stored needed to do a Hopf point continuation. This includes the
initializers and the Hopf point curve definition file.
• LimitCycle
Here are all files stored needed to do a limit cycle continuation. This includes the
initializers and the limitcycle curve definition file.
• LimitCycleCodim2
This directory contains the routines that compute the normal form coefficients of codim2
bifurcations of limit cycles. Computing these coefficients reliably requires that compu-
tations are done to high precision and symbolic derivatives up to fifth order are used.
• LimitPoint
Here are all files stored needed to do a limitpoint continuation. This includes the
initializers and the limitpoint curve definition file.
• LimitPointCycle
Here are all files stored needed to do a fold bifurcation of limit cycles continuation. This
includes the initializers and limitpoint of cycles curve definition files.
• MultiLinearForms
Here are files stored needed to compute high-order derivatives of systems and normal-
form coefficients of bifurcations.
• NeimarkSacker
Here are all files stored needed to do a torus bifurcation of limit cycles continuation.
This includes the initializers and torus curve definition files.
• PeriodDoubling
Here are all files stored needed to do a period doubling bifurcation continuation. This
includes the initializers and perioddoubling curve definition files.
• SBML
This directory contains the routines that allow to import SBML (Systems Biology
Markup Language) systems into MatCont. These routines require an update and/or
revision.
• Systems
Here all system definition files and files related to computed curves are stored, except
those odefiles which are hard-stored in the directory Testruns/TestSystems because
these are used in the testruns. To avoid confusion, it is not recommended to use the
names of the odefiles that are present in Testruns/TestSystems.
30
• Testruns
Here are some example testruns stored which can be executed from the matlab com-
mand line. They can be used to run the examples described in this manual and to test
if everything is working correctly. The necessary odefiles are stored in the subdirectory
TestSystems.
We note that 12 directories are dedicated to a specific curve type, namely:

Equilibrium, LimitCycle, LimitPoint, Hopf,
LimitPointCycle, PeriodDoubling, NeimarkSacker, BranchPoint,
BranchPointCycle, Homoclinic, Heteroclinic and HomoclinicSaddleNode.
Three other directories, namely HomotopyHet, HomotopySaddle and HomotopySaddleNode
are initialization directories to a specific curve type.
The remaining 8 directories are not dedicated to a specific curve type.
The only files which are not in any of these directories are init.m, cpl.m, matcont.fig
and matcont.m. The function init must be called before using the command-line toolbox
so that matlab can find all the needed functions. cpl is used to plot the results of the
continuation in CL MatCont (for details see §3.9). matcont.m is the start-up file of the
GUI version MatCont, and matcont.fig is the related figure-file.
3.8 Global Structures

In general the user does not need to know much about the use of global structures in Mat-
Cont but advanced users can find it useful since a lot of additional information can be hidden
in these structures. MatCont uses a structure cds (continuation descriptor structure) which
is global in the continuer routine cont.m and in all (12) curve definition files and it carries
information back and forth between them (cf. Figure 4). Also, cds has a field cds.options
with the fields described in §3.4. It contains default values for these options but overrides
them with the values in the options-structure that is passed with the call of cont.m. Other
fields of cds.options are filled via the interaction with the curve definition file which is also
passed with the call to cont.m
Next, MatCont uses more specific global structures, namely:
• eds (equilibrium descriptor structure): this structure is global in the curve definition
file equilibrium.m and in the initializers for the continuation of equilibria.
• lpds (limit point descriptor structure): this structure is global in the curve definition
file limitpoint.m and in the initializers for the continuation of limit points.
• hds (Hopf descriptor structure): this structure is global in the curve definition file
hopf.m and in the initializers for the continuation of Hopf points.
• bpds (branch point descriptor structure): this structure is global in the curve definition
file branchpoint.m and in the initializers for the continuation of branch points.
• lds (limitcycle descriptor structure): this structure is global in the curve definition files
limitcycle.m, limitpointcycle.m, perioddoubling.m, neimarksacker.m, and also
in branchpointcycle.m, in the initializers for the continuation of limit cycles, their
codimension 1 bifurcations and branch points.
31
• homds (homoclinc descriptor structure): this structure is global in the curve definition
files homoclinic.m, homoclinicsaddlenode.m and in the initializers for the continua-
tion of orbits to saddle or to saddle-node.
• hetds (heteroclinc descriptor structure): this structure is global in the curve definition
file heteroclinic.m and in the initializers for the continuation of heteroclinic orbits.
These specific structures carry information that is collected or computed in the initializers
to the curve definition files in which they are global. In the initializers themselves other
specific structures can also be global, depending on the data that are used in the initializer.
For example, in the initializer init H LC not only lds and cds are global but also eds and
hds. Indeed, Hopf points can be detected on equilibrium curves or be taken from Hopf curves.
3.9 Plots in MatCont

When MatCont is handled from the command line, then of course all MATLAB plot func-
tionalities can be used. However, MatCont provides two specific plot functions, namely
cpl.m and plotcycle.m.
cpl(x, v, s, e) makes a two or three dimensional plot. This routine automatically
places labels at the singularity points.The first three arguments must be the output of a
continuation run. The fourth argument e is optional. e is an array whose elements define
which coordinates of the system are used. Therefore e must have either 2 components (2D-
plot) or 3 components (3D-plot). If e is not given and x has 2 (respectively 3) components,
then a 2D-plot (3D-plot,respectively) is drawn. In all other cases an error message will be
generated.
The use of plotcycle is discussed in §7.3.
32
4 The odefile of a dynamical system
4.1 Structure and construction of an odefile
A solution curve must be initialized before doing a continuation. Each curve file has its own
initializers and in the case of dynamical systems the initializers typically use an odefile where
the ode is defined. An odefile (short for system definition file) contains at least the following
sections:
init, fun eval, jacobian, jacobianp, hessians, hessiansp, der3, der4 , der5.
An odefile may also contain one or more sections that describe user functions.
There is a variety of options to create such odefiles. In all cases users should avoid using
loops, conditional statements, or other specific programming constructions in the system
definition.
First, an odefile can be defined by simply using the Matlab editor (or, in fact, any text
editor). This is likely to lead to errors and therefore not encouraged if the matlab symbolic
toolbox is available and symbolic derivatives are desirable.
A second option is to use the GUI version of MatCont by defining the problem in
the ’System’ window and then choosing the options “symbolically” or “numerically” to use
symbolic or numeric (= finite difference) derivatives, respectively. The user functions can also
be added in the GUI of MatCont, using the menu ‘User function’.
From version 6.7 on MatCont provides another possibility to create odefiles. It is in fact
a shortcut to using the GUI version of MatCont. As an example we will study the Rössler
chaotic system: 
 ẋ = −y − z
ẏ = x + Ay

ż = Bx − Cz + xz,
where (x, y, z) are the phase variables, and (A, B, C) are the parameters.
An odefile can be created by calling SysGUI.new.
This opens a System window, which contains several fields and buttons. To identify the
system, type for example
Roessl
in the Name field (it must be one word).

Input names of the Coordinates: X,Y,Z, and the Parameters: AA,BB,CC.
If shown, select symbolic generation of the 1st order derivatives by pressing the corre-
sponding radio-button 1 .
Finally, in the large input field, type the RHS of the truncated normal form map as
X’=-Y-Z
Y’=X+AA*Y
Z’=BB*X-CC*Z+X*Z
1
If the MATLAB Symbolic Toolbox is present, there will be buttons indicated ’symbolically’. The first-
order derivatives are used in some of the integration algorithms, the first- and second-order derivatives are
used in the continuation, while the third-order derivatives are employed in the normal form computations.
The derivatives of fourth and fifth order are only used in the normal form computations of some codimension
2 bifurcations.
33
Figure 5: Specifying a new model.
34
Avoid typical mistakes:
• Make sure the multiplication is written explicitly with ∗.
• Specify the right hand sides in the same order as the coordinates.
It is best not to add comma’s or semicolons after the equations. Now the System window
should look like in Figure 5, and you can press the OK button. Two new files will be created in
the Systems directory of MatCont, namely the odefile Roessl.m and a mat-file Roessl.mat.
The odefile can be edited later on by calling SysGUI.edit(@name) where name is the name
of an existing odefile. User functions can be added by calling SysGUI.userfunctions(@name).
See Figure 6.
Figure 6: Adding a user function.
We note that if state variables, parameters or user functions are added or deleted then
this constitutes another dynamical system. So either all computed data should be deleted or
ignored, or the name of the system should be changed. The last option is recommended, in
particular when the GUI is used.
A special point on a bifurcation curve that is specified by a user function has a structure
as follows:
35
s.index index of the detected singular point defined by the user function.
s.label a string that is in UserInfo.label, label of the singularity.
s.data an empty tangent vector, values of the test and user functions in
the singular point.
s.msg a string that is set in UserInfo.name.
We now give an odefile for the Rössler system using symbolic derivatives of all orders up
to three and including the userfunction x − 1. Note that we have returned to state variables
denoted by x, y, z and parameters A, B, C.
function out = Roessl

out{1} = @init;
out{2} = @fun_eval;
out{3} = @jacobian;
out{4} = @jacobianp;
out{5} = @hessians;
out{6} = @hessiansp;
out{7} = @der3;
out{8} = [];
out{9} = [];
out{10}= @xis1;
% --------------------------------------------------------------------------
function dydt = fun_eval(t,kmrgd,A,B,C)
dydt=[-kmrgd(2)-kmrgd(3);
kmrgd(1)+A*kmrgd(2);
B*kmrgd(1)-C*kmrgd(3)+kmrgd(1)*kmrgd(3);];
% --------------------------------------------------------------------------
function [tspan,y0,options] = init
y0=[0,0,0];
options = odeset(’Jacobian’,handles(3),’JacobianP’,handles(4),...
’Hessians’,handles(5),’HessiansP’,handles(6));
handles = feval(Roessl);
tspan = [0 10];
% --------------------------------------------------------------------------
function jac = jacobian(t,kmrgd,A,B,C)
jac=[ 0 , -1 , -1 ; 1 , A , 0 ; B + kmrgd(3) , 0 , kmrgd(1) - C ];
% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,A,B,C)
jacp=[ 0 , 0 , 0 ; kmrgd(2) , 0 , 0 ; 0 , kmrgd(1) , -kmrgd(3) ];
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,A,B,C)
hess1=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 1 ];
hess2=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
hess3=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 1 , 0 , 0 ];
hess(:,:,1) =hess1;
36
hess(:,:,2) =hess2;
hess(:,:,3) =hess3;
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,A,B,C)
hessp1=[ 0 , 0 , 0 ; 0 , 1 , 0 ; 0 , 0 , 0 ];
hessp2=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 1 , 0 , 0 ];
hessp3=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , -1 ];
hessp(:,:,1) =hessp1;
%---------------------------------------------------------------------------
function tens3 = der3(t,kmrgd,A,B,C)
tens31=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens32=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens33=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens34=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens35=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens36=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens37=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens38=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens39=[ 0 , 0 , 0 ; 0 , 0 , 0 ; 0 , 0 , 0 ];
tens3(:,:,1,1) =tens31;
tens3(:,:,1,2) =tens32;
tens3(:,:,1,3) =tens33;
tens3(:,:,2,1) =tens34;
tens3(:,:,2,2) =tens35;
tens3(:,:,2,3) =tens36;
tens3(:,:,3,1) =tens37;
tens3(:,:,3,2) =tens38;
tens3(:,:,3,3) =tens39;
%---------------------------------------------------------------------------
%---------------------------------------------------------------------------
function userfun1=xis1(t,kmrgd,A,B,C)
userfun1=kmrgd(1)-1;
We observe the following:
• The state variables are collected in a vector called kmrgd.
• Internally the names of the parameters are extended so that, e.g., “C” is replaced by
“par C”. This is done to avoid clashes with the symbolic toolbox in which certain
names are protected. This does not affect the use of the odefile since the parameters
are not referred to by name. If the user writes his own odefile then this precaution is
not necessary.
37
• The init function is called when performing a time integration in the GUI mode of
MatCont. It is of no use in the command line version. See §5.1.3 for a command line
example.
4.2 Handling auxiliary functions in the construction of the odefile

If complicated expressions arise several times in the definition of a dynamical system then
is often useful to give them a name. A simple example is the following system (a catalytic
oscillator, (60):

 ẋ = 2q1 z 2 − 2q5 x2 − q3 xy
ẏ = q2 z − q6 y − q3 xy (31)

ṡ = q4 z − kq4 s
where z = 1 − x − y − s. When SysGUI.new is used to create the odefile of the system then
is not necessary to “expand” the auxiliary function for z. The system can be introduced as
z=1-x-y-s
x’=2*q1*z^2-2*q5*x^2-q3*x*y
y’=q2*z-q6*y-q3*x*y
z’=q4*z-k*q4*s
We note that the definition of the auxiliary functions must precede the definiton of the
right hand sides.
The corresponding odefile is cataloscill.m in the directory Testruns/TestSystems.
4.3 Access to function and Jacobian values

Once an odefile is created, it can be used to compute the values of the odefunction and its
derivatives with respect to state variables and parameters on the command line. As an ex-
ample, consider the Morris-Lecar system MyML in the directory Testruns/TestSystems. This
system has two state variables and two parameters. We compute the odefunction values and
derivatives with respect to the state vaariables and parameters for the state vector [0; 0] and
parameters (30,10) by executing the script testodefile, available in the directory Testruns:
MyMLSystem=MyML
MyMLfunction=MyMLSystem{2}
MyMLjacobian=MyMLSystem{3}
MyMLjacobianp=MyMLSystem{4}
MyMLfunction(0,[0;0],30,10)
MyMLjacobian(0,[0;0],30,10)
MyMLjacobianp(0,[0;0],30,10)
The output is as follows:
MyMLSystem =
1 x 9 cell array
38
{@init} {@fun_eval} {@jacobian} {@jacobianp} {@hessians}
{@hessiansp} {0x0 double} {0x0 double} {0x0 double}
MyMLfunction =
function_handle with value:
@fun_eval
MyMLjacobian =
@jacobian
MyMLjacobianp =
@jacobianp
ans =
33.1953
0.0167
ans =
1.8282 -128.0000
0.0013 -0.0694
ans =
0.2000 0
0 -0.0013
4.4 Time integration using the odefile

Time integration (= simulation) of the dynamical system can be done on the command line
as in the following example:
OPTIONS=[];
39
hls = MyML;
[t,y] = ode45(hls{2},[0 1000],[0 0],OPTIONS,30,10);
x0 = y(end,:)’;
In this example the system MyML is integrated over the interval [0, 1000] starting with the
input vector [0, 0] and with parameters 30, 10 in that order. The integration is performed
with the integration routine ode45 with the defaults (OPTIONS=[]) in Matlab. The output
vector is placed in x0. More details are given in Chapter 5.
40
5 Time integration and Poincaré maps
matlab provides a suite of ODE-solvers. To create a combined integration-continuation
environment, we have made them all accessible in MatCont and added two new ones, ode78
and ode87. ode78 is an explicit Runge-Kutta method that uses 7th order Fehlberg formulas
[Fehlberg 1969]. This is a 7-8th-order accurate integrator, therefore the local error normally
expected is O(h9 ). ode87 is a Runge-Kutta method that uses 8-7th order Dorman and Prince
formulas. See [Prince and Dorman 1981]. This is a 8th-order accurate integrator therefore
the local error normally expected is O(h9 ). In MatCont the choice of the solver is made
via the Integrator window that is opened automatically when the curve type O (Orbit) is
chosen. We note that the Starter window has a SelectCycle button that allows to start the
continuation of periodic orbits from curves computed by time integration.
5.1 Time integration

MatCont uses the matlab representation of ODEs, which can also be used by the matlab
ODE solvers. The user of Cl MatCont is therefore also able to use his/her model within
the standard ODE solver without rewriting the code. In MatCont, also backward time inte-
gration is possible using the standard ODE solvers. To make them accessible in MatCont,
some output functions and properties were created.
5.1.1 Solver output properties

The solver output properties allow one to control the output that the solvers generate. Mat-
Cont detects whether extra output is needed by looking at the number of open output
windows (2D, 3D or numeric). If extra output is required, MatCont sets the OutputFcn
property to the output function, integplot which is passed to an ODE solver by options
= odeset(’OutputFcn’, @integplot), otherwise it is set to the function integ prs. The
output function must be of the form status = integplot(t, y, flag, p1 , p2 ,...). The
solver calls this function after every successful integration step with the following flags. Note
that the syntax of the call differs with the flag. The function must respond appropriately:
• init: The solver calls integplot(tspan,y0,’init’) before beginning the integration,

to allow the output function to initialize. This part initializes the output windows
and does some initializations to speed up the further processing of the output. It also
launches the window that makes it possible to interactively ”stop/pause/resume” the
computations.
• within MatCont it is possible to define the number of points (npoints) after which
output is needed. Because the solver calls status = integplot (t,y) after each in-
tegration step, the number of calls to integplot and npoints does not correspond.
Therefore this part is divided into two parts. t contains points where output was gener-
ated during the step, and y is the numerical solution at the points in t. If t is a vector,
the i-th column of y corresponds to the i-th element of t. The output is produced ac-
cording to the Refine option. integplot must return a status output value of 0 or 1. If
status = 1, the solver halts integration. This part also handles the ”stop/pause/resume”
interactions.
41
• done: The solver calls integplot([],[],’done’) when integration is completed to
allow the output function to perform any cleanup chores. The stop-window is also
deleted.
Setting the OutputFcn property to the output function, integ prs, reduces the output time.
It only allows to interactively pause, resume and stop the integration.
5.1.2 Jacobian matrices

Stiff ODE solvers often execute faster if the Jacobian matrix, i.e. the matrix of partial deriva-
tives of the RHS that defines the differential equations, is provided. The Jacobian matrix
pertains only to those solvers for stiff problems (ode15s, ode23s, ode23t, ode23tb), for
which it can be critical for reliability and efficiency. If the user does not provide a code to
calculate the Jacobian matrix, these solvers approximate it numerically using finite differ-
ences. Supplying an analytical Jacobian matrix often increases the speed and reliability of
the solution for stiff problems. If the Jacobian matrix is provided (symbolically) in the odefile,
MatCont stores as property the function handle of the Jacobian, otherwise this handle is
set empty.
The standard matlab odeget and odeset only support Jacobian matrices of the RHS
w.r.t. the state variables. However, we do need derivatives with respect to the parameters
for the continuation. To compute normal form coefficients, it is also useful to have higher-
order symbolic derivatives. To overcome this problem, MatCont contains new versions of
odeget and odeset, which support Jacobian matrices with respect to parameters and higher-
order partial derivatives w.r.t state variables. The new routines are compatible with the ones
provided by matlab.
5.1.3 The Steinmetz-Larter example

In the peroxidase-oxidase reaction a peroxidase catalyzes an aerobic oxidation:
2Y H2 + O2 + 2H + → 2Y H + + 2H2 O,
where Y H2 (N ADH) is a general electron donor. Under the proper conditions this reaction
yields oscillations in the concentrations of Y H2 , O2 and various reaction intermediates. Stein-
metz and Larter [Steinmetz and Larter 1991] derived the following system of four coupled
nonlinear differential rate equations:


 Ȧ = −k1 ABX − k3 ABY + k7 − k−7 A,

Ḃ = −k1 ABX − k3 ABY + k8 ,
(32)

 Ẋ = k1 ABX − 2k2 X 2 + 2k3 ABY − k4 X + k6 ,

Ẏ = −k3 ABY + 2k2 X 2 − k5 Y,
where A, B, X, Y are state variables and k1 , k2 , k3 , k4 , k5 , k6 , k7 , k8 , and k−7 are parameters.

Although highly simplified, the model is able to reproduce the three modes of simple, chaotic,
and bursting oscillations found in experiments. State and parameter values of a point on
a NS cycle in (32) are given in Table 4. The normal form coefficient of the NS cycle is
−1.406017e − 006. Since it is negative, we expect stable tori nearby.
To start a time integration from this NS point, with a slightly supercritical parameter
value, namely k7 = 0.7167, we can use the commands
42
Variable Value Parameter Value Parameter Value
A 1.8609653 k1 0.1631021 k5 1.104
B 25.678306 k2 1250 k6 0.001
X 0.010838258 k3 0.046875 k7 0.71643356
Y 0.094707061 k4 20 k−7 0.1175
k8 0.5
Table 4: State and parameter values of a point on an NS- cycle in the Steinmetz-Larter model.
X0=[1.8609653;25.678306;0.010838258;0.094707061];
P0=[0.1631021;1250;0.046875;20;1.104;0.001;0.71643356;0.1175;0.5];
P0(7)=0.7167;
hls=feval(@SteLar);
options=odeset(’RelTol’,1e-8);
[t,y]=ode45(hls{2},[0,3000],X0,options,0.1631021,1250,0.046875,20,1.104,
0.001,0.7167,0.1175,0.5);
size(t)
size(y)
plot(t(100000:339000),y(100000:339000,2));
’press any key’
pause
plot(y(100000:339000,3),y(100000:339000,2));
After a transient the time-series exhibits modulated oscillations with two frequencies near
the original limit cycle (see Fig. 7). This is a motion on a stable two-dimensional torus that
arises from the Neimark-Sacker bifurcation. The commands can be executed by running the
testrun testtimeinteg.m. The total number of computed points is 339161.
26.4
26.2
26
B
25.8
25.6
25.4
25.2
500 1000 1500 2000 2500 3000
t
(a) Modulated oscillations. (b) Orbits on a stable 2-torus.
Figure 7: Movement on a stable torus in a 4-dimensional space as a function of time in the

B− direction and as projected on the (B, X) phase plane.
The testrun testtimeintegJacobian.m
X0=[1.8609653;25.678306;0.010838258;0.094707061];
P0=[0.1631021;1250;0.046875;20;1.104;0.001;0.71643356;0.1175;0.5];
P0(7)=0.7167;
hls=feval(@SteLar);
43
options=odeset(’RelTol’,1e-8);
options = odeset(’Jacobian’,hls(3),’JacobianP’,hls(4),
’Hessians’,hls(5),’HessiansP’,hls(6));
[t,y]=ode45(hls{2},[0,3000],X0,options,0.1631021,1250,0.046875,20,
1.104,0.001,0.7167,0.1175,0.5);
size(t)
size(y)
plot(t(100000:325000),y(100000:325000,2));
’press any key’
pause
plot(y(100000:325000,3),y(100000:325000,2));
is a small modification of testtimeinteg.m in which the solver has access to the Jacobian
matrix. The graphical output is similar, but the number of computed points is now 328197.
5.2 Poincaré section and Poincaré map

A Poincaré section is a surface in phase space that cuts across the flow of a dynamical
system. It is a carefully chosen (in general, curved) surface in the phase space that is crossed
by almost all orbits. It is a tool developed by Poincaré for visualization of the flow in more
than two dimensions. The Poincaré section has one dimension less than the phase space
and the Poincaré map transforms the Poincaré section onto itself by relating two consecutive
intersection points, say uk and uk+1 . We note that only those intersection points count,
which come from the same side of the section. The Poincaré map is invertible because one
gets un from un+1 by following the orbit backwards. A Poincaré map turns a continuous-time
dynamical system into a discrete-time one. If the Poincaré section is carefully chosen no
information is lost concerning the qualitative behaviour of the dynamics. For example, if the
system is being attracted to a limit cycle, one observes dots converging to a fixed point in the
Poincaré section.
5.2.1 Poincaré maps in Cl MatCont

When computing an orbit (t, y(t)) in Matlab an event can be defined as going through a
zero of a given scalar function G(t, y). If G does not explicitly depend on time in an au-
tonomous dynamical system, this feature can be used to detect Poincaré intersections. One
does this by setting the Events property to a function handle, e.g. @events, creating a
function [value,isterminal,direction] = events(t,y) and calling
[t,Y,TE,YE,IE] = solver(odefun,tspan,y0,options).
For the i-th event function:
• value(i) is the value of the function.
• isterminal(i) = 1 if the integration is to terminate at a zero of this event function

and 0 otherwise.
• direction(i) = 0 if all zeros are to be computed (the default), +1 if only the zeros
are needed where the event function increases, and -1 if only the zeros where the event
function decreases.
44
Corresponding entries in TE, YE, and IE return, respectively, the time at which an event oc-
curs, the solution at the time of the event, and the index i of the event function that vanishes.
Example
Let a set of two event functions be introduced by defining the function testEV:
function [value,isterminal,direction]= testEV(t,y,varargin)

value=[y(1)-0.2;y(2)-0.3];
isterminal=zeros(2,1);
direction=ones(2,1);
end
This event function requires that the system is at least two-dimensional and defines two
events, namely y(1) = 0.2 and y(2) = 0.3. The integration will not be terminated if an
event is detected and all zeros will be detected regardless of the direction of y (increasing or
decreasing).
We now consider the system adaptx in the directory Testruns/TestSystems of MatCont
with three state variables and two parameters. By runnning the script testPoincare:
TSTP=@testEV;
OPTIONS = odeset(’RelTol’,1e-8,’Events’,TSTP);
hls = adaptx;
[t,y,TE,YE,IE] = ode45(hls{2},[0 300],[0.3 0.5 -0.1],OPTIONS,1,0.8);
x0 = y(end,:);
[t,y,TE,YE,IE] = ode45(hls{2},[0 10],x0,OPTIONS,1,0.8);
we integrate the system adaptx from 0 to 300 and then further from 300 to 310 starting from
the point [0.3; 0.5; −0.1] with parameter values 1, 0.8 and we detect the points in the second
run where y(1) = 0.2 or y(2) = 0.3. The output is given by
TE =
5.6334
6.6979
YE =
-0.2247 0.3000 0.3131

0.2000 0.4213 -0.1059
IE =
2
1
The test is run by typing testPoincare in the command line.
45
6 Equilibrium continuation
This example will show how to continue an equilibrium of a differential equation defined in
a standard matlab ODE file. Furthermore, this example illustrates the detection, location,
and processing of singularities along an equilibrium curve.
6.1 Mathematical definition

Consider a differential equation
du
= f (u, α), u ∈ Rn , α ∈ R f : Rn+1 → Rn . (33)
dt
We are interested in an equilibrium curve, i.e. f (u, α) = 0. The defining function is therefore:
F (x) = f (u, α) = 0 (34)
with x = (u, α) ∈ Rn+1 . Denote by v ∈ Rn+1 the tangent vector to the equilibrium curve at
x.
6.2 Initialization by time integration

If not much previous knowledge is available on a given dynamical system then a common way
to start the study is to compute orbits and see whether some of them converge to (stable)
equilibria or periodic orbits. An example for the case of equilibria is given in the testrun
testmymlPRC.m.
6.3 Bifurcations and their normal form coefficients

In continuous-time systems there are two generic codim 1 bifurcations that can be detected
along the equilibrium curve (no derivations will be done here; for more detailed information
see [25]):
• fold, also known as limit point. We will denote this bifurcation by LP
• Hopf-point, denoted by H
The equilibrium curve can also have branch points. These are denoted with BP. To detect
these singularities, we first define 3 test functions:

Fx
φ1 (u, α) = det , (35)
vT
  
0
 (2fu (u, α) ⊙ In ) w1  ... 
φ2 (u, α) =   T \
 0 
 , (36)
v1 d
1 n+1
φ3 (u, α) = vn+1 , (37)
where ⊙ is the bialternate matrix product and v1 , w1 are n(n−1) 2 vectors chosen so that the
square matrix in (36) is non-singular. (the bialternate matrix product was introduced by C.
Stéphanos [36]; see [22] for details). Using these test functions we can define the singularities:
46
• BP: φ1 = 0
• H: φ2 = 0
• LP: φ3 = 0, φ1 6= 0
A proof that these test functions correctly detect the mentioned singularities can be found in
[25]. Here we only notice that φ2 = 0 not only at Hopf points but also at neutral saddles, i.e.
points where fx has two real eigenvalues with sum zero. So, the singularity matrix is:
 
0 − −
S= − 0 −  (38)
1 − 0
For each detected limit point, the corresponding quadratic normal form coefficient is com-
puted:
1
a = pT fuu [q, q], (39)
2
where fu q = fuT p = 0, q T q = 1, pT q = 1. Mathematically, the limit point is nondegenerate
(i.e. the equilibrium branch looks locally like a parabola) if and only if a 6= 0. In practice,
because of round-off errors the computed a is always nonzero. So it is more important to
check if the testfunction φ3 changes sign. (In a cusp point (CP) a = 0 but this will not be
detected on a branch of equilibria since it is a codimension 2 phenomenon, see §8.1)
At a Hopf bifurcation point, the first Lyapunov coefficient is computed by the formula
1
l1 = Re pT fuuu [q, q, q̄] − 2fuu [q, Fu−1 fuu [q, q̄]] + fuu [q̄, (2iωIn − fu )−1 fuu [q, q]] , (40)
2
where fu q = iωq, q T q = 1, fuT p = −iωp, p̄T q = 1.
The first Lyapunov coefficient is quite important. If l1 < 0 then the Hopf bifurcation
is supercritical, i.e. within the center manifold on one side of the bifurcation only stable
equilibria exist, on the other side unstable equilibria coexist with stable periodic orbits. If
l1 > 0 then the Hopf bifurcation is subcritical, i.e. within the center manifold on one side
of the bifurcation stable equilibria coexist with unstable periodic orbits, on the other side
only unstable equilibria exist. (In a Generalized Hopf point (GH) l1 = 0 but this will not be
detected on a branch of equilibria since it is a codimension 2 phenomenon, see §8.2)
6.3.1 Branch point locator

The location of Hopf and limit points usually does not cause problems. However, the location
of branch points can give problems. The region of attraction of the Newton type continuation
method which is used, has the shape of a cone (see [2]). In the localisation process we cannot
assume to stay in this cone. This difficulty can be avoided by introducing p ∈ Rn and β ∈ R
and considering the extended system:


 f (u, α) + βp = 0

fuT (u, α)p = 0
(41)

 pT fα (u, α) = 0

pT p − 1 = 0
47
We solve this system by Newton’s method with initial data β = 0 and p : fuT p = µp where
µ is the real eigenvalue with smallest norm. A branch point (u, α) corresponds to a regular
solution (u, α, 0, p) of system (41) (see [3],p. 165). We note that the second order partial
derivatives (Hessian) of f with respect to u and α are required.
The tangent vector at the singularity is also computed here. This is related to the pro-
cessing of the branch point (computing the direction of the secondary branch).
6.4 Equilibrium initialization

Naively, one would start the continuation immediately with:
[x,v,s,h,f]=cont(@equilibrium, x0, v0, opt)
However, the equilibrium curve file has to know :
• which ode file to use,
• the values of all state variables,
• the values of all parameters,
• which parameter is active.
All this information can be supplied by one of the following two starting functions. Both
functions return an initial point x0 as well as its tangent vector v0.
• [x0,v0]=init EP EP(@odefile, x, p, ap)
This routine stores its information in a global structure eds (see also Figure 4). The
result of init EP EP is a vector x0 with the state variables and the active parameter
and a vector v0 that is empty. Here odefile is the ode-file (system-definition file) to be
used, x is a vector containing the values of the state variables. p is the vector containing
the current values of the parameters and ap is the active parameter. The full listing
of the equilibrium initializer can be found in the file ’Equilibrium/init EP EP.m’ of the
toolbox.
• [x0,v0]=init BP EP(@odefile, x, p, s, h)
Calculates an initial point for starting a new branch from a branch point detected on
an equilibrium curve. This routine stores its information in a global structure eds (see
also Figure 4). Here odefile is the ode-file (system-definition file) to be used, x is a
vector containing the values of the state variables returned by a previous equilibrium
curve continuation. p is the vector containing the current values of the parameters
and h contains the value of the initial amplitude. The full listing of the equilibrium
initializer and the curve definition can be found in the files ’Equilibrium/init BP EP.m’
and ’equilibrium.m’ of the toolbox, respectively.
6.5 Bratu example

The first example we will look at is a 4-point discretization of the Bratu-Gelfand BVP [22].
This model is defined as follows:
x′ = y − 2x + aex (42)
′ y
y = x − 2y + ae (43)
48
The system is specified as
bratu.m
1 function out = bratu
2 out{1} = @init;
3 out{2} = @fun_eval;
4 out{3} = @jacobian;
5 out{4} = @jacobianp;
6 out{5} = @hessians;
7 out{6} = @hessiansp;
8 out{7} = [];
9 out{8} = [];
10 out{9} = [];
11 out{10}= @userf1;
12
13 end
14
15 % --------------------------------------------------------------------------
16 function dydt = fun_eval(t,kmrgd,a)
17 dydt = [ -2*kmrgd(1)+kmrgd(2)+a*exp(kmrgd(1));
18 kmrgd(1)-2*kmrgd(2)+a*exp(kmrgd(2)) ];
19
20 % --------------------------------------------------------------------------
21 function [tspan,y0,options] = init
22 tspan = [0; 10];
23 y0 = [0;0];handles = feval(@bratu)
24 options = odeset(’Jacobian’,handles(3),’JacobianP’, ’handles(4)’,...
25 ...,’Hessians’,handles(5), ’Hessiansp’,handles(6));
26 % --------------------------------------------------------------------------
27 function jac = jacobian(t,kmrgd,a)
28 jac = [ -2+a*exp(kmrgd(1)) 1
29 1 -2+a*exp(kmrgd(2)) ];
30
31 % --------------------------------------------------------------------------
32 function jacp = jacobianp(t,kmrgd,a)
33
34 jacp = [ exp(kmrgd(1))
35 exp(kmrgd(2)) ];
36
37 % --------------------------------------------------------------------------
38 function hess = hessians(t,kmrgd,a)
39 hess1=[[a*exp(kmrgd(1)),0];[0,0]];
40 hess2=[[0,0];[0,a*exp(kmrgd(2))]];
41 hess(:,:,1) = hess1;
42 hess(:,:,2) = hess2;
43
44 % --------------------------------------------------------------------------
45 function hessp = hessiansp(t,kmrgd,a)
46 hessp1=[[exp(kmrgd(1)),0];[0,exp(kmrgd(2))]];
47 hessp(:,:,1) = hessp1;
48
49 %---------------------------------------------------------------------------
60 function userfun1 = userf1(t,kmrgd,a)
61 userfun1 = a-0.2;
62
bratu.m
49
As seen above, a user function is defined to detect all points where a = 0.2. This system
has an equilibrium at (x, y, a) = (0, 0, 0) which we will continue with respect to a. We first
compute 50 points and then extend the curve with another 50 points.
global cds
p=[0];ap=[1];
[x0,v0]=init_EP_EP(@bratu,[0;0],p,ap);
opt=contset;
opt=contset(opt,’MaxNumPoints’,50);
opt=contset(opt,’Singularities’,1);
opt=contset(opt,’Userfunctions’,1);
UserInfo.name=’userf1’;
UserInfo.state=1;
UserInfo.label=’u1’;
opt=contset(opt,’UserfunctionsInfo’,UserInfo);
[x,v,s,h,f]=cont(@equilibrium,x0,[],opt);
[x,v,s,h,f]=cont(x,v,s,h,f,cds);
cpl(x,v,s,[3 1 2]);
The above computations can be done by running testbratu.m in the directory Testruns.
The output in the command window is as follows:
>> testbratu
first point found
tangent vector to first point found
label = u1, x = ( 0.259174 0.259174 0.200002 )
label = LP, x = ( 1.000001 1.000001 0.367879 )
a=3.535537e-01
Neutral saddle
label = H , x = ( 2.000000 2.000000 0.270671 )
label = u1, x = ( 2.542639 2.542639 0.200000 )
elapsed time = 0.4 secs
npoints curve = 50
start computing extended curve
label = BP, x = ( 3.000000 3.000000 0.149361 )
npoints curve = 100
We note that in the first continuation two zeros of the user function (label u1) were detected,
as well as a limit point (label LP) and a neutral equilibrium (label H). A neutral equilibrium
is an equilibrium with two real eigenvalues with sum zero.
In the extension of the run at (x, y, a) ≈ (3.0; 3.0; 0.15) the system has a branch point
(label BP).
cpl(x,v,s,[3 1 2]) plots a 3D-plot with the parameter a on the x-axis and the first and
second state variable on the y- and z-axes, respectively. The labels of the plot are changed
manually by the following commands:
50
Figure 8: Equilibrium curve of bratu.m with a branch point
• press ‘Insert’ → ‘XLabel’ →, ‘a’
• press ‘Insert’ → ‘YLabel’ →, ‘x’
• press ‘Insert’ → ‘ZLabel’ →, ‘y’
The resulting curve is plotted in Figure 8.
To select the branch point the output s is used. Since the first and last points are also
treated as singular, the array of structures s has 7 components and the data concerning
the branch point are at s(6). To switch to another branch at the detected branch point, we
select that branch point and we use the starter init BP EP. We start a forward and backward
continuation from this point:
testbratu;
x1=x(1:2,s(6).index);
p(ap)=x(3,s(6).index);
[x0,v0]=init_BP_EP(@bratu,x1,p,s(6),0.01);
opt=contset(opt,’InitStepsize’,0.0001);
[x1,v1,s1,h1,f1]=cont(@equilibrium,x0,v0,opt);
cpl(x1,v1,s1,[3 1 2]);
opt=contset(opt,’Backward’,1);
[x2,v2,s2,h2,f2]=cont(@equilibrium,x0,v0,opt);
cpl(x2,v2,s2,[3 1 2]);
Note that p is a vector containing the initial values of all parameters. The above compu-
tations can be done by running testbratu2.m in the directory Testruns. The output in the
command window is as follows:
>> testbratu2
first point found
51
Figure 9: Equilibrium curve of bratu.m with new branches rooted in the branch point
label = u1, x = ( 0.259174 0.259174 0.200002 )

label = LP, x = ( 1.000001 1.000001 0.367879 )
a=3.535537e-01
label = H , x = ( 2.000000 2.000000 0.270671 )
label = u1, x = ( 2.542639 2.542639 0.200000 )

npoints curve = 50
start computing extended curve
label = BP, x = ( 3.000000 3.000000 0.149361 )

npoints curve = 100
first point found

npoints curve = 50
first point found
label = BP, x = ( 3.000000 3.000000 0.149361 )

npoints curve = 50
We note that the branch point is ”discovered” a second time during the backward contin-
uation of the secondary branch that is rooted at the branch point. All computed curves are
plotted together in Figure 9.
52
7 Continuation of limit cycles
7.1 Mathematical definition
Consider the following differential equation
du
= f (u, α) (44)
dt
with u ∈ Rn and α ∈ R. A periodic solution with period T satisfies the following system
du
dt = f (u, α) . (45)
u(0) = u(T )
For simplicity the period T is treated as a parameter resulting in the system

du
dτ = T f (u, α) . (46)
u(0) = u(1)
If u(τ ) is its solution then the shifted solution u(τ + s) is also a solution to (46) for any value
of s. To select one solution, a phase condition is added to the system. The complete BVP
(boundary value problem) is
 du
 dτ − T f (u, α) = 0
u(0) − u(1) = 0 (47)
 R1
0 hu(t), u̇ old (t)idt = 0
where u̇old is the derivative of a previous solution. A limit cycle is a closed phase orbit
corresponding to this periodic solution.
7.2 Discretization of a limit cycle

In MatCont limit cycles are discretized using orthogonal collocation [6], the same way as
it was done in AUTO [12]; the left hand side of the resulting system is the defining function
F (u, T, α) for limit cycles.
In practice this means that the normalized time interval [0, 1] is divided in a number ntst
(number of test intervals) of intervals with variable lengths; the (ntst+1) endpoints of these
intervals form the coarse mesh. Each interval is further subdivided in a number ncol (number
of collocation points) of subintervals with equal lengths. So altogether there are (ntst × ncol
+ 1) fine mesh points. The state variable values are collected in the points of the fine mesh.
The limit cycle itself is then approximated by a continuous piecewise polynomial which is a
polynomial of degree ncol in each of the ntst coarse mesh intervals.
If the number of phase variables is denoted by nphase then altogether
(ntst × ncol × nphase)
state variable values are stored, since by periodicity the values in the endpoints 0 and 1 must
be identical. But the number of continuation variables is
((ntst × ncol + 1 ) × nphase)+2
since the state variable values in both 0 and 1 are stored and the period T and a free sys-
tem parameter must be included. They are stored in the output vector x of a limit cycle
continuation in that order, after all state variable values.
53
The coarse mesh is adapted after each number of Adapt continuation points, cf §3.4.9. We
note that ntst and ncol are input arguments of the routines that initialize the continuation of
limit cycles but are not explicitly found in the output of the continuation; however, they are
preserved as fields in the global structure lds.
7.3 Plotting the output of a continuation of limit cycles

A continuation run of limit cycles has the standard form
[x,v,s,h,f]=cont(@limitcycle,x0,v0,opt);
Here x is a matrix in which each column corresponds to a computed orbit. The last
element of each column contains the value of the active parameter. Cl MatCont provides
only limited possibilities to visualise the output. Namely, it provides the routine plotcycle.m
which is to be called in the form
plotcycle(x,v,s,e)
where e is either a two- or a three-dimensional row vector.
• If e=[e1,e2] one must have 1 ≤ e1, e2 ≤ nphase and a 2D plot will be generated in
which the projections of the orbits on the (e1, e2) coordinate plane are shown. The first
and last curves, and the curves with singulaties are in red.
• If e=[e1,e2,e3] then e1 must be the index of the active parameter, i.e. e1=size(x,1)
and one must have 1 ≤ e2, e3 ≤ nphase. Then a 3D plot will be generated in which the
projections of the orbits on the (e2, e3) coordinate plane are shown versus the active
parameter. The first and last curves, and the curves with singulaties are in red.
7.4 Initialization by time integration

The initialization of a continuation of limit cycles is a nontrivial issue. Probably the most
often used method is to start from a Hopf point that is detected on a curve of equilibria, see
§7.6. Another powerful method is to compute orbits and see whether some of them converge
to (stable) periodic orbits. If such an orbit is found (the easiest way is by graphical inspection)
then one uses an (approximation to) a point of the orbit and integrates again over a time
interval that is somewhat larger than the period of the orbit but smaller that twice the period.
The routine initOrbLC.m in the directory LimitCycle then can initialize the continuation
of periodic orbits with a user-chosen free parameter, with the period of orbit as the second
free parameter. An example is provided in the testrun testselectcycle.m in the directory
Testruns:
OPTIONS = [];
hls = adaptx;
OPTIONS=odeset(’RelTol’,1e-8);
[t,y] = ode45(hls{2},[0 300],[0.3 0.5 -0.1],OPTIONS,1,0.8);
x1 = y(end,:);
[t,y] = ode45(hls{2},[0 10],x1,OPTIONS,1,0.8);
figure
54
0.5
0.4
0.3
0.2
0.1
-0.1
-0.2
-0.3
-0.4
-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
Figure 10: A closed curve in the adaptx - system.
plot(y(:,1),y(:,2))
p=[1;0.8];
ap=[2];
tolerance=1e-2;
[x0,v0]=initOrbLC(@adaptx,t,y,p,ap,20,4,tolerance);
opt=contset;
%opt=contset(opt,’TSearchOrder’,0);
%opt=contset(opt,’Backward’,1);
[xlcc,vlcc,slcc,hlcc,flcc]=cont(@limitcycle,x0,v0,opt);
figure
axes
plotcycle(xlcc,vlcc,slcc,[size(xlcc,1) 1 2]);
In this script adaptx is the dynamical system (54) in §7.7 whose odefile adaptx.m is stored
in the directory Testruns/TestSystems. It has three state variables and two parameters.
Starting from the point (0.3; 0.5; −0.1) and with parameter values (1, 0.8) it is integrated over
the time span 300. The endpoint is called x1 and from this point another time integration
over the shorter span 10 is performed. It is then checked graphically that this new integration
contains a closed circle, see Figure 10 for a plot in a two-dimensional space.
The call
[x0,v0]=initOrbLC(@adaptx,t,y,p,ap,20,4,tolerance);
now initializes the continuation of limit cycles. Here t is a column vector whose entries are
the time points (between 0 and 10) and y is a matrix whose rows contain the coordinates of
the points computed along the orbit. Also, p = (1; 0.8) is the parameter vector and ap = [2]
contains the indices of the free parameters, so in this case the second parameter is free. Also,
the number of test functions is 20 and the number of collocation points in the discretization
55
0.4
0.2
y
0
-0.2
0.8
0.6 0.8
0.4 0.75
0.2 0.7
0 0.65
-0.2 0.6
x -0.4 0.55 alpha
Figure 11: 50 limit cycles in the adaptx - system.
of the limit cycle is 4. Finally, the tolerance is a threshold for accepting an initial part of the
computed orbit as an approximation to a limit cycle.
The last command plots the continuation of a branch of limit cycles with 50 computed
limit cycles in the state space, see Figure 11 (the axis labels were added manually).
Note: when the lines
%opt=contset(opt,’TSearchOrder’,0);
in testselectcycle are made active, then the initialization of the continuation of limit
cycles involves a decreasing order of the unit vectors in continuation space instead of an
increasing order. This accidentally leads to a continuation in the opposite direction; therefore
cds.options.Backward is set to reverse the continuation.
7.5 Bifurcations of limit cycles

On a limit cycle curve the following bifurcations can occur
• Branch Point of Cycles, this will be denoted as BPC
• Period Doubling, denoted as PD
• Fold, also known as Limit Point of Cycles, this will be denoted as LPC
• Neimark-Sacker, this will be denoted as NS
The test function for the Period Doubling bifurcation is defined by the following system

 v̇(τ ) − T fu (u, α)v(τ ) + Gϕ(τ ) = 0
v(0) + v(1) = 0 (48)
 R1
0 hψ(τ ), v(τ )idτ = 1
here ϕ and ψ are so-called bordering vector-functions [25], see [15] for details on the im-
plementation. The system is discretized using orthogonal collocation and solved using the
56
standard matlab sparse system solver. The solution component G ∈ R of this system is the
test function and equals zero when there is a Period Doubling bifurcation.
The Fold bifurcation is detected in the same way as the Fold bifurcation of equilibria, the
last component of the tangent vector (the α component) is used as the test function.
The Neimark-Sacker bifurcation is detected by monitoring the eigenvalues of the mon-
odromy matrix for the cycle. The monodromy matrix is computed like in auto by a block
elimination in the discretized form of the Jacobian of (47).
BPC cycles are not generic in families of limit cycles, but they are common in the case of
symmetries, if the branch parameter is also the continuation parameter. CL MatCont uses
a strategy that requires only the solution of linear systems; it is based on the fact that in a
symmetry-breaking BPC cycle MD has rank defect two, where MD is the square matrix MD ,
obtained from the discretized form of the Jacobian of (47). To be precise, if h ∈ C 1 ([0, 1], IRn ),
then
ḣ − T fx (x(t), α)h
Mh = ,
h(0) − h(1)
and
(ḣ − T fx (x(t), α)h)dc
MD (h)dm = ,
h(0) − h(1)
where ()dm and ()dc denote discretization in mesh points and in collocation points, respectively.
Therefore we border MD with two additional rows and columns to obtain
 
MD w 1 w 2
MDbb =  v1∗ 0 0 ,
∗
v2 0 0
so that MDbb is nonsingular in the BPC cycle. Then we solve the systems
   
ψ11 ψ12 0(N m+1)n 0(N m+1)n
MDbb  gBP C11 gBP C12  =  1 0 ,
gBP C21 gBP C22 0 1
where ψ11 , ψ12 have (N m + 1)n components, and gBP C11 , gBP C12 , gBP C21 , and gBP C22 , are
scalar test functions for the BPC. In the BPC cycle they all vanish.
The singularity matrix is
 
0 0 0 0 − − − −
 − − − − − 0 − − 
S=  − − − − − − 0 − 
 (49)
− − − − 1 − 1 0
The first row corresponds to the BPC. It contains 4 zeros which indicates that gBP C11 , gBP C12 ,
gBP C21 , and gBP C22 should vanish. The last row corresponds to the NS. Because we have to
exclude that all four testfunctions of the BPC are zeros, we introduce an extra testfunction
which corresponds to the norm of these four testfunctions. A NS is detected if this norm is
nonzero,the testfunction for the fold is nonzero and the testfunction for the NS is equal to
zero.
57
7.5.1 Branch Point Locator
The location of BPC points in the non-generic situation (i.e. where some symmetry is present)
as zeros of the test functions is numerically suspect because no local quadratic convergence
can be guaranteed. This difficulty can be avoided by introducing an additional unknown
β ∈ R and considering the minimally extended system:
 dx

 dt − T f (x, α) + βp1 =0
 x(0) − x(1) + βp2 =0
R1 (50)

 hx(t), ẋold (t)idt + βp3 = 0
 0
G[x, T, α] =0
where G is defined as in (84) and [pT1 pT2 p3 ]T is the bordering vector [w01 ; w02 ; w03 ]T in (86).
We solve this system with respect to x, T, α and β by Newton’s method with initial β = 0. A
branch point (x, T, α) corresponds to a regular solution (x, T, α, 0) of system (50) (see [3],p.
165). We note, however that the second order partial derivatives (Hessian) of f with respect
to x and α are required. The tangent vector v1st at the BPC singularity is approximated as
v1st = v1 +v
2
2
where v1 is the tangent vector in the continuation point previous to the BPC
and v2 is the one in the next point.
7.5.2 Normal form coefficients

The numerical methods used in MatCont to compute normal form coefficients of codimen-
sion 1 bifurcations of limit cycles are discussed in [28].
The periodic normal form at the Limit Point of Cycles (LPC) bifurcation is
( dτ
dt = 1 − ξ + aξ 2 + · · · ,
dξ
(51)
2
dt = bξ + · · · ,
where τ ∈ [0, T ], ξ is a real coordinate on the center manifold that is transverse to the limit-
cycle, a, b ∈ R, and dots denote nonautonomous T -periodic O(ξ 3 )-terms. For each detected
LPC point the normal form coefficient b is computed. If b 6= 0 then the LPC bifurcation is
nondegenerate, i.e., the cycle manifold is quadratically tangential to the hyperplane orthogo-
nal to the parameter direction. In practice, by round-off errors b is never zero, so it is better
to check if the testfunction changes sign. (for b = 0 the LPC degenerates to a cusp of cycles
(CP) but this will not be detected on a branch of limit cycles since it is a codimension 2
phenomenon.)
The periodic normal form at the Period Doubling (PD) bifurcation is
( dτ
dt = 1 + aξ 2 + · · · ,
dξ
(52)
3
dt = cξ + · · · ,
where τ ∈ [0, 2T ], ξ is a real coordinate on the center manifold that is transverse to the limit
cycle, a, c ∈ R, and dots denote nonautonomous 2T -periodic O(ξ 4 )-terms. The coefficient c
determines the stability of the period doubled cycle in the center manifold and is computed
during the processing of each PD point.
If c < 0 then the PD bifurcation is supercritical, i.e., within the center manifold on one
side of the bifurcation only stable cycles exist and on the other side unstable cycles coexist
58
with stable period-doubled cycles. If c > 0 the PD bifurcation is subcritical, i.e., within
the center manifold on one side of the bifurcation stable cycles coexist with unstable period-
doubled cycles and on the other side only unstable cycles exist. (for c = 0 the PD degenerates
to Generalized Period Doubling bifurcation (GPD) but this will not be detected on a branch
of limit cycles since this is a codimension 2 phenomenon).
The periodic normal form at the Neimark-Sacker (NS) bifurcation is
( dτ
dt = 1 + a|ξ|2 + · · · ,
dξ
(53)
iθ 2
dt = T ξ + dξ|ξ| + · · · ,
where τ ∈ [0, T ], ξ is a complex coordinate on the center manifold that is complementary

to τ , a ∈ R, d ∈ C, and dots denote nonautonomous T -periodic O(|ξ|4 )-terms. The critical
coefficient d in the periodic normal form for the NS bifurcation is computed during the
processing of a NS point.
If Re d < 0 then the bifurcation is supercritical, i.e. in the center manifold the limit cycle
is stable on one side of the bifurcation point and unstable on the other side and the unstable
cycles coexist with stable invariant tori. If Re d > 0 then the bifurcation is subcritical, i.e. in
the center manifold the limit cycle is stable on one side of the bifurcation point and unstable
on the other side and the stable cycles coexist with unstable invariant tori. (for Re d = 0 the
NS bifurcation degenerates to a Chenciner bifurcation (CH) but this will not be detected on
a branch of limit cycles since it is a codimension 2 phenomenon).
7.6 Limitcycle initialization

For limit cycles the same problems occur as with equilibria, a limit cycle continuation can’t
be done by just calling the continuer as
[x,v,s,h,f]=cont(@limitcycle, x0, v0, opt)
The limit cycle curve file has to know
• which parameter is active,
• the number of mesh and collocation points to use for the discretization.
Also an initial cycle x0 has to be known. All this information can be supplied using any of the
following three starting functions. All three return an initial cycle x0 as well as its tangent
vector v0.
• [x0,v0]=init H LC(@odefile, x, p, ap, h, ntst, ncol)

Calculates an initial cycle from a Hopf point detected on an equilibrium curve. Here
odefile is the ode-file to be used. x is a vector containing the values of the state vari-
ables returned by a previous equilibrium curve continuation. p is the vector containing
the current values of the parameters. ap is the active parameter, h contains the value of
the initial amplitude and ntst and ncol are the number of mesh and collocation points
to be used for the discretization.
59
• [x0,v0]=init BPC LC(@odefile, x, v, s, ap, ntst, ncol,h)
Calculates an initial cycle for starting a secondary cycle from a BPC detected on a pre-
vious calculated limit cycle. odefile, ap, ntst and ncol are the same as for init H LC.
x, v and s are here the x, v and s as returned by a previous limit cycle continuation
and h contains the value of the initial amplitude.
• [x0,v0]=init LC LC(@odefile, x, v, s, par, ap, ntst, ncol)

This starter can be used to start a limit cycle continuation from a previous limit cycle
continuation. odefile, ap, ntst and ncol are the same as for init H LC. x and v are
here the x and v as returned by a previous limit cycle continuation. s is the special point
structure for the point from where to start the new continuation. par is the parameter
vector that is to be used for the new continuation. If this is the same as in the cycle
returned by the continuer (arguably the natural situation), then this is the same as
s.data.parametervalues
• [x0,v0]=init PD LC(@odefile, x, s, ntst, ncol,h)

This starter calculates an initial double period cycle from a period doubling bifurcation
point. Here h is a distance parameter which can be set to zero. We note also that there
is no active parameter ap input argument.
7.7 Adaptive control example

For this example the following system from adaptive control was used in a feedback control
system, as described in [20], [21] and further used in [25] (Example 5.4, p. 178):

 ẋ = y
ẏ = z (54)

ż = −αz − βy − x + x2
This system is introduced in the testruns of MatCont under the name adaptx. It has a
Hopf point at the origin for α = 1 and β = 1 as detected in the testrun testadapt.m.
From this Hopf point an initial cycle is calculated using the starter init H LC. The results
of the continuation are plotted using the plot function plotcycle(x,v,s,e) (see Figure 12).
This function plots the cycles x. e is an array whith either 2 or 3 elements for 2-dimensional
and 3-dimensional plotting respectively. Its entries must be indices of state variables or active
parameters in x. The index of the active parameter is size(x,1)
>> init;
>> [x0,v0]=init_EP_EP(@adaptx,[0;0;0],[-10;1],[1]);
>> opt = contset; opt = contset(opt,’Singularities’,1);
>> [x,v,s,h,f]=cont(@equilibrium,x0,[],opt);
first point found
label = H , x = ( 0.000000 0.000000 0.000000 1.000002 )
First Lyapunov coefficient = -3.000001e-001

npoints curve = 300
>> x1=x(1:3,s(2).index);p=[x(end,s(2).index);1];
60
>> [x0,v0]=init_H_LC(@adaptx,x1,p,[1],1e-6,20,4);
>> opt = contset(opt,’MaxNumPoints’,200);
>> opt = contset(opt,’Multipliers’,1);
>> opt = contset(opt,’Adapt’,1);
>> [xlc,vlc,slc,hlc,flc]=cont(@limitcycle,x0,v0,opt);
first point found
Limit point cycle (period = 6.283185e+000, parameter = 1.000000e+000)
Normal form coefficient = -1.306379e+000
Branch Point cycle(period = 6.283185e+000, parameter = 9.999996e-001)
Period Doubling (period = 6.364071e+000, parameter = 6.303020e-001)
Normal form coefficient = -4.267675e-002
Neimark-Sacker (period = 6.433818e+000, parameter = 1.895460e-008)
Neutral saddle
Period Doubling (period = 6.364071e+000, parameter = -6.303020e-001)
Normal form coefficient = 4.268472e-002

npoints curve = 200
>> plotcycle(xlc,vlc,slc,[size(xlc,1) 1 2]);
We note that xlc is a 245 × 200 matrix; each column corresponds to a computed limit
cycle and gives the coordinates of all points of the fine mesh, i.e. 243 = (20 × 4 + 1) × 3 values,
plus the period T as the 244 − th component and the value of the active parameter α as the
245 − th component.
The x-axis contains the active parameter, the y-axis the first state variable x and the
z-axis the second state variable y. This run can be tested by the statement testadapt1 (the
axis labels have to be set manually). If you run only this example, do not forget to execute
init init statement first.
We note that the Limit point cycle and Branch point cycle detected in this run are
degenerate: they reduce to the Hopf point itself. The Neimark-Sacker bifurcation is, in
reality, a Neutral Saddle, i.e. it is an unstable periodic orbit with two real multipliers whose
product is 1.
From the first Period Doubling bifurcation detected a limit cycle continuation of the nearby
double period cycle is started. First, an initial cycle and its tangent vector are calculated
using the starter init PD LC. The continuation is done using the standard continuer and the
result is plotted using the plotcycle function (see figure 13).
>> [x1,v1]=init_PD_LC(@adaptx,xlc,slc(4),40,4,1e-6);
>> opt=contset(opt,’MaxNumPoints’,250);
>> [xlc2,vlc2,slc2,hlc2,flc2]=cont(@limitcycle,x1,v1,opt);
first point found
61
Figure 12: Computed limit cycle curve
Neutral saddle
Limit point cycle (period = 1.103168e+001, parameter = -4.494912e-002)
Normal form coefficient = 1.311327e+002
Neimark-Sacker (period = 1.076785e+001, parameter = -1.076152e-009)
Neutral saddle
Limit point cycle (period = 1.103169e+001, parameter = 4.494912e-002)
Neutral saddle

npoints curve = 250
>> plotcycle(xlc2,vlc2,slc2,[size(xlc2,1) 1 2]);
This run can be tested by the statement testadapt2.
62
Figure 13: Computed limit cycle curve started from a Period Doubling bifurcation
7.8 The phase response curve

The phase response curve of a limit cycle, or PRC, is a curve, defined over the period of the
cycle, that expresses, at each time of that period, the effect of a small input vector on the
cycle. In experimental circumstances, this may correspond to injected current, to the addition
of more chemical agents, etc. A positive value means that the current cycle is shortened in
time, a negative value means that the period is prolonged.
The PRC, as it is generally computed, is exact for infinitesimally small input vectors. In
practice the maximum norm of the input vector would depend on the needed accuracy and
the values of the system’s state variables.
The derivative phase response curve or dPRC also has some very important applications.
For the concrete use of PRC and dPRC in synchronization studies in neural modeling, we
refer to [23].
Cl MatCont supports the computation of the PRC and dPRC of limit cycles during
continuation, using the method described in [24]. The standard method, which uses numerical
integration of the adjoint system, was implemented in XPPAUT [17].
The use in MatCont is easy: before starting the actual limit cycle continuation, the user
can specify whether he wants to compute the PRC, dPRC or both, and he needs to indicate
the input vector used. When a scalar is given as input, then the vector has this scalar as first
entry and all other entries are zero. Then in separate plotting windows, for each computed
step in limit cycle continuation, the PRC and/or dPRC are computed and plotted. The
computed values are saved in the output f-array of the continuation.
63
To illustrate the use in CL MatCont, we here supply code that does a LC continuation
experiment in the Morris Lecar system (the odefile is called MyML.m and is located in the
directory Testruns/TestSystems), and also computes the PRC and dPRC of the system:
clear
global x v s h f opt
OPTIONS=[];
hls = MyML;
[t,y] = ode45(hls{2},[0 1000],[0 0],OPTIONS,30,10);
x0 = y(end,:)’;
opt=contset;opt=contset(opt,’Singularities’,1);
opt = contset(opt,’MaxStepsize’,10);
opt=contset(opt,’MaxNumpoints’,2500);
[x1,v1]=init_EP_EP(@MyML,x0,[30;6],[1]);
x1=x(1:2,s(5).index);p=[x(end,s(5).index);6];
[x0,v0]=init_H_LC(@MyML,x1,p,[1],1e-6,40,4);
opt = contset(opt,’Multipliers’,0);
opt = contset(opt,’Adapt’,1);
opt = contset(opt,’MaxStepsize’,5);
opt = contset(opt,’FunTolerance’,1e-6);
opt = contset(opt,’VarTolerance’,1e-6);
opt = contset(opt,’PRC’,1);
opt = contset(opt,’dPRC’,1);
opt = contset(opt,’Input’,1);
opt = contset(opt,’MaxNumPoints’,100);
[xlc2,vlc2,slc2,hlc2,flc2]=cont(@limitcycle,x0,v0,opt);
fvector=flc2(:,100);
PRC10=fvector(42:202);
dPRC10=fvector(203:363);
plot(PRC10,’r’);
hold on
plot(dPRC10,’b’);
We note that the first 41 components (number of test intervals plus one) of the f vector
contain the coarse mesh points of the discretization. The next 161 components (number of
fine mesh points, i.e. number of test intervals times number of collocation points plus 1)
contain the values of the PRC in the fine mesh points. The next 161 components contain the
values of the dPRC in the fine mesh points.
This script is in the MatCont directory Testruns under the name testmymlPRC. The
output is:
>> testmymlPRC
first point found
64
0.015
0.01
0.005
PRC and dPRC
-0.005
-0.01
-0.015
-0.02
-0.025
0 20 40 60 80 100 120 140 160 180
Phase
Figure 14: Phase response curve in the Morris-Lecar system: red=PRC, blue=dPRC
label = H , x = ( -28.700773 0.018189 43.312018 )

First Lyapunov coefficient = 6.461827e-03
label = LP, x = ( -26.127744 0.024296 43.740592 )
a=9.306236e-03
label = LP, x = ( -10.804133 0.126584 34.546930 )
a=-2.105189e-03
label = H , x = ( 7.947868 0.555741 197.796289 )
First Lyapunov coefficient = 8.882425e-04
Current step size too small (point 1134)
npoints curve = 1134
first point found

npoints curve = 100
The code also produces Figure 14.

To further illustrate the computation of the PRC and the dPRC in CL MatCont, we
here supply code that does a LC continuation experiment in the adaptx system (the odefile is
called adaptx.m and is located in the directory Testruns/TestSystems), and also computes
the PRC and dPRC of the system:
[x0,v0]=init_EP_EP(@adaptx,[0;0;0],[-10;1],[1]);
opt=contset;opt=contset(opt,’Singularities’,1);
65
x1=x(1:3,s(2).index);p=[x(end,s(2).index);1];
[x0,v0]=init_H_LC(@adaptx,x1,p,[1],1e-6,20,4);
opt = contset(opt,’Multipliers’,1);
opt = contset(opt,’Adapt’,1);
[xlc,vlc,slc,hlc,flc]=cont(@limitcycle,x0,v0,opt);
par=slc(end).data.parametervalues;
[x1,v1] = init_LC_LC(@adaptx, xlc, vlc, slc(end), par, 1, 20, 4);
opt = contset(opt,’PRC’,1);
opt = contset(opt,’dPRC’,1);
opt = contset(opt,’Input’,1);
[xlc1,vlc1,slc1,hlc1,flc1]=cont(@limitcycle,x1,v1,opt);
fvector=flc1(:,20);
plot(fvector(22:102),’r’)
hold on
plot(fvector(103:183),’b’)
We note that the first 21 components (number of test intervals plus one) of the f vector
contain the coarse mesh points of the discretization. The next 81 components (number of fine
mesh points, i.e. number of test intervals times number of collocation points plus 1) contain
the values of the PRC in the fine mesh points. The next 81 components contain the values
of the dPRC in the fine mesh points.
This script is in the MatCont directory Testruns under the name testadaptPRC. The
output is:
>> testadaptPRC
first point found
label = H , x = ( 0.000000 0.000000 0.000000 1.000002 )

npoints curve = 300
first point found

npoints curve = 10
first point found
66
1.5
0.5
PRC and dPRC
-0.5
-1
-1.5
0 10 20 30 40 50 60 70 80 90
Phase
Figure 15: Phase response curve in the adaptx system: red=PRC, blue=dPRC

npoints curve = 20
The code also produces Figure 15.
67
8 Continuation of codim 1 bifurcations
MatCont computes branches of codim 1 bifurcations of equilibria (namely, fold LP and
Hopf H) and of limit cycles (fold of cycles LPC, period doubling PD and Neimark-Sacker
NS). On these branches codimension 2 bifurcations of equilibria (CP,BT,ZH,HH,GH) and
of limit cycles (CPC,R1,R3,R4,CH,LPNS,PDNS,R2, NSNS,LPPD and GPD) are detected
and located, cf. the graph of adjacency in Figure 1.1. For the numerical methods used in
the computation of the normal form coefficients, but also for the meaning and use of these
coefficients we refer to [25] in the equilibrium case and to [9] in the limit cycle case. We note
that the MatCont directory LimitCycleCodim2 is exclusively devoted to the computation
of normal form coefficients of codim 2 bifurcations of limit cycles.
8.1 Fold Continuation

8.1.1 Mathematical definition
In the toolbox fold curves are computed by minimally extended defining systems cf. [22],
§4.1.2. The fold curve is defined by the following system

f (u, α) = 0,
(55)
g(u, α) = 0,
where (u, α) ∈ Rn+2 , while g is obtained by solving

fu (u, α) wbor v 0n
T = , (56)
vbor 0 g 1
and wbor , vbor ∈ Rn are chosen such that the matrix in (56) is nonsingular. An advantage of
this method is that the derivatives of g can be obtained easily from the derivatives of fu (u, α):
gz = −wT (fu )z v
where z is a state variable or an active parameter and w is obtained by solving

T
fu (u, α) vbor w 0n
T = , (57)
wbor 0 g 1
This method is implemented in the curve definition file limitpoint.m.
8.1.2 Bifurcations along a fold curve

In continuous-time systems there are four generic codim 2 bifurcations that can be detected
along a fold curve:
• Bogdanov - Takens. We will denote this bifurcation by BT
• Zero - Hopf point, denoted by ZH
• Cusp point, denoted by CP
• Branch point, denoted by BP
68
To detect these singularities, we first define bp + 3 test functions, where bp is the number of
branch parameters:
• φ1 = wT v (cf. formula (39))
•   
0
 (2fu (u, α) ⊙ In ) w1  ... 
φ2 (u, α) = 
 \ 
 0 
 (58)
v1T d
1 n+1
• φ3 = wT fuu [v, v]
• φi = wT fβi (u, α)
In these expressions v, w are the vectors computed in (56) and(57) respectively, ⊙ is the
bialternate matrix product, v1 , w1 are n(n−1)
2 vectors chosen so that the square matrix in (58)
is non-singular, and βi (branch parameters) are components of α. The singularity matrix for
bp = 0 is:  
0 0 −
S= 1 0 −  (59)
− − 0
The number of branch parameters is not fixed. If the number of branch parameters is 2
then this matrix has two more rows and columns. This singularity matrix is automatically
extended:  
0 0 − − −
 1 0 − − − 
 
S=  − − 0 − − 

 − − − 0 − 
− − − − 0
8.1.3 Fold initialization

The most natural way to start a fold curve continuation in the toolbox is from a limit point.
As in the case of equilibria, a fold curve continuation cannot be done by just calling the
continuer as:
[x,v,s,h,f]=cont(@limitpoint, x0, v0, opt).
The fold curve file has to know
• which odefile to use,
• which parameters are active,
• the values of all parameters.
To initialize the continuation one first gives the following command:

[x0,v0]=init LP LP(@odefile, xnew, p, ap (, bp)optional ).
In this command xnew must be a vector that contains the values of the state variables.
p must contain the current values of all the parameters and ap must be the indices of
69
the 2 active parameters. In the most natural situation where x is the matrix returned
by the previous equilibrium curve continuation one starts to build xnew by the command
xnew=x(1:nphase,s(i).index), where nphase is the number of state variables and s(i) is
the special point structure of the detected fold point on the equilibrium curve continuation.
Next, the command p(ap old)=x(end,s(i).index); replaces the old value of the free pa-
rameter in the previous run by the newly found parameter p. odefile specifies the ode-file to
be used. bp are the optional indices of the branch parameters. It also works without entering
a value for this field.
MatCont provides four other initializers which allow to continue a fold curve from a
codim 2 equilibrium bifurcation, namely init BP LP.m, init BT LP.m, init CP LP.m and
init ZH LP.m. In fact, these initializers are only added for ease of use: they refer back to
init LP LP.m.
8.1.4 Adaptation
It is possible to adapt the problem while generating the fold curve. This call updates the
auxiliary variables used in the defining system of the computed branch. The bordering vectors
vbor and wbor may require updating since they must at least be such that the matrices in (56)
and (57) are nonsingular. Updating is done by replacing vbor and wbor by the normalized
vectors v, w computed in (56) and (57) respectively.
8.1.5 Example: a catalytic oscillator

For this example the following system (a catalytic oscillator) is used

 ẋ = 2q1 z 2 − 2q5 x2 − q3 xy
ẏ = q2 z − q6 y − q3 xy (60)

ṡ = q4 z − kq4 s
where z = 1 − x − y − s and the parameters q1 , q2 , q3 , q4 , q5 , q6 , k are introduced in that order

in the odefile cataloscill in the MatCont directory Testruns/TestSystems.
The starting vector x0 and its tangent vector v0 are calculated from the following equi-
librium curve continuation (q2 is free).
p=[2.5;2.204678;10;0.0675;1;0.1;0.4];
ap1=[2];
[x0,v0]=init_EP_EP(@cataloscill,[0.001137;0.891483;0.062345],p,ap1);
opt=contset;
opt=contset(opt,’MaxStepSize’,0.025);
cpl(x,v,s,[4,1]);
This set of commands is run by executing the file testequilcataloscill in the Mat-
Cont directory Testruns. The command line output is the following:
>> testequilcataloscill
first point found
70
0.16
0.14
0.12
0.1
0.08 H
x
0.06 LP
0.04
LP
0.02 H
0
0.8 1 1.2 1.4 1.6 1.8 2 2.2
q2
Figure 16: Computed equilibrium curve in the catalytic oscillar system (60)

label = H , x = ( 0.016357 0.523973 0.328336 1.051558 )
First Lyapunov coefficient = 1.070259e+01
label = LP, x = ( 0.024717 0.450257 0.375018 1.042049 )
a=-1.166509e-01
label = LP, x = ( 0.054030 0.302241 0.459807 1.052200 )
a=1.346534e-01
label = H , x = ( 0.077929 0.233063 0.492149 1.040991 )

npoints curve = 78
The results are plotted using the plot function cpl where the fourth argument selects the
fourth and first components of the solution which are the parameter q2 and the coordinate x.
The results can be seen in Figure 16 where the axes labels are added manually.
We now add a forward and a backward fold continuation from the first LP detected on
the previous equilibrium curve; q2 and k are free in both runs.
p=[2.5;2.204678;10;0.0675;1;0.1;0.4];
ap1=[2];
opt=contset;
cpl(x,v,s,[4,1]);
’press any key’
71
pause
x1=x(1:3,s(3).index);
p(ap1)=x(end,s(3).index);
[x0,v0]=init_LP_LP(@cataloscill,x1,p,[2,7]);
[x2,v2,s2,h2,f2]=cont(@limitpoint,x0,v0,opt);
hold on;
cpl(x2,v2,s2,[4,1]);
’press any key’
pause
hold on
cpl(x3,v3,s3,[4,1]);
This set of commands is run by executing the file testLPcataloscill in the MatCont
directory Testruns. The command line output is the following:
>> testLPcataloscill
first point found
label = H , x = ( 0.016357 0.523973 0.328336 1.051558 )
label = LP, x = ( 0.024717 0.450257 0.375018 1.042049 )
a=-1.166509e-01
label = LP, x = ( 0.054030 0.302241 0.459807 1.052200 )
a=1.346534e-01
label = H , x = ( 0.077929 0.233063 0.492149 1.040991 )

npoints curve = 78
ans =
’press any key’
first point found

label = CP , x = ( 0.035941 0.352004 0.451371 1.006408 0.355991 )
c=3.627908e-01
label = BT , x = ( 0.115909 0.315467 0.288437 1.417627 0.971397 )
(a,b)=(8.378438e-02, 2.136279e+00)

npoints curve = 78
ans =
72
0.16
0.14
0.12 BT
0.1
0.08 H
x
0.06 LP
0.04 CP
LP
0.02 H BT
0
0.8 1 1.2 1.4 1.6 1.8 2 2.2
q2
Figure 17: Computed equilibrium and fold curves in the catalytic oscillator system (60)
’press any key’
first point found

label = BT , x = ( 0.016337 0.638410 0.200456 1.161199 0.722340 )
(a,b)=(-4.822577e-02, -1.937637e+00)

npoints curve = 78
The results are plotted in Figure 17, where, as usual, the axis labels are added manually.
8.2 Hopf Continuation

In the MatCont / CL MatCont toolbox Hopf curves are computed by minimally extended
defining systems, cf. [22] §4.3.4. The Hopf curve is defined by the following system

 f (u, α) = 0,
gi j (u, α, k) = 0 (61)
 11
gi2 j2 (u, α, k) = 0

g11 g12
with the unknowns u, α, k, (i1 , j1 , i2 , j2 ) ∈ {1, 2} and where g = is obtained by
g21 g22
solving
fu2 + kIn Wbor V 0n,2
T = , (62)
Vbor O G I2
73
where fu has eigenvalues ±iω, ω > 0, k = ω 2 and Vbor , Wbor ∈ Rn×2 are chosen such that the
matrix in (62) is nonsingular. i1 , j1 , i2 , j2 , Vbor and Wbor are auxiliary variables that can be
adapted. This method is implemented in the curve definition file hopf.m.
8.2.2 Bifurcations along a Hopf curve

In continuous-time systems there are four generic codim 2 bifurcations that can be detected
along a Hopf curve:
• Bogdanov - Takens. We will denote this bifurcation by BT
• Zero - Hopf point, denoted by ZH
• Double - Hopf point, denoted by HH
• Generalized Hopf point, denoted by GH
To detect these singularities, we first define 4 test functions:
• φ1 = k
• φ2 = det(fu )
  
0
 (2fu (u, α) ⊙ In ) w1  ... 
• φ3 = 
 \
 0


v1T d
1 n+1
• φ4 = l1 (first Lyapunov coefficient, see (40)),

n(n−1)
where v1 , w1 are carefully constructed and updated 2 × 2 matrices.
In this case the singularity matrix is:
 
0 0 − −
 1 0 − − 
S=  − − 0
 (63)
− 
− − − 0
8.2.3 Hopf initialization

The most natural way to start the continuation of a Hopf bifurcation curve is to start it
from a Hopf point, typically found on an equilibrium curve. The continuation of Hopf points
cannot simply be started by using the following command:
[x,v,s,h,f]=cont(@hopf, x0, v0, opt)
The Hopf curve file has to know
• the values of all state variables,
• which parameters are active,
• the values of all parameters.
74
Also an initial point x0 has to be known. All this information can be supplied using the
following command:
[x0,v0]=init H H(@odefile, xnew, p, ap).
The input arguments odefile, xnew, p, ap are built exactly as in init LP LP. The output
vector x0 consists of xnew extended with the free parameters and paramter k from (62).
MatCont provides four other initializers which allow to continue a Hopf curve by starting
from a codim 2 equilibrium point. These are init BT H.m, init GH H.m, init HH H.m and
init ZH H.m These initializers are added merely for ease of use since they refer back to
init H H.m. However, we note that two Hopf curves pass through a Hopf point. So it can
be necessary to specify the tangent vector to the desired Hopf curve.
8.2.4 Adaptation
It is possible to adapt the problem while generating the Hopf curve. This call updates
the auxiliary variables used in the defining system of the computed branch. The bordering
matrices Vbor and Wbor may require updating since they must at least be such that the matrix
in (62) is nonsingular. Updating of Vbor and Wbor is done exactly as in the LP case. Updating
of i1 , j1 , i2 , j2 is done in such a way that the linearized system of (61) is as well-conditioned
as possible.
8.2.5 Example
For this example we again use the system (60). The starting vector is calculated from the
first Hopf bifurcation point detected in the equilibrium continuation example (see 8.1.5). q2
and k are free.
p=[2.5;2.204678;10;0.0675;1;0.1;0.4];
ap1=[2];
opt=contset;
cpl(x,v,s,[4,1]);
’press any key’
pause
x1=x(1:3,s(3).index);
[x0,v0]=init_LP_LP(@cataloscill,x1,p,[2,7]);
hold on;
cpl(x2,v2,s2,[4,1]);
’press any key’
pause
75
hold on
cpl(x3,v3,s3,[4,1]);
’press any key’
pause
x1=x(1:3,s(2).index);
[x0,v0]=init_H_H(@cataloscill,x1,p,[2 7]);
opt=contset;
opt=contset(opt,’MaxStepsize’,0.1);
opt=contset(opt,’backward’,0);
[x4,v4,s4,h4,f4]=cont(@hopf,x0,v0,opt);
hold on;
cpl(x4,v4,s4,[4 1]);
This set of commands is executed by the file testLPHopfcataloscill.m in the directory
Testruns. The command line output is
>> testLPHopfcataloscill
first point found
label = H , x = ( 0.016357 0.523973 0.328336 1.051558 )
label = LP, x = ( 0.024717 0.450257 0.375018 1.042049 )
a=-1.166509e-01
label = LP, x = ( 0.054030 0.302241 0.459807 1.052200 )
a=1.346534e-01
label = H , x = ( 0.077929 0.233063 0.492149 1.040991 )

npoints curve = 78
ans =
’press any key’
first point found

label = CP , x = ( 0.035941 0.352004 0.451371 1.006408 0.355991 )
c=3.627908e-01
label = BT , x = ( 0.115909 0.315467 0.288437 1.417627 0.971397 )
(a,b)=(8.378438e-02, 2.136279e+00)

npoints curve = 78
76
ans =
’press any key’
first point found

label = BT , x = ( 0.016337 0.638410 0.200456 1.161199 0.722340 )
(a,b)=(-4.822577e-02, -1.937637e+00)

npoints curve = 78
ans =
’press any key’
first point found

label = GH, x = ( 0.018022 0.368238 0.497968 0.891319 0.232487 0.003324 )
l2=-7.768956e+02
label = GH, x = ( 0.064311 0.211095 0.554870 0.924255 0.305879 0.003512 )
l2=-2.401240e+02
label = BT, x = ( 0.115909 0.315467 0.288437 1.417627 0.971396 0.000000 )
(a,b)=(8.378437e-02, 2.136280e+00)
label = BT, x = ( 0.016337 0.638410 0.200456 1.161199 0.722339 0.000000 )
(a,b)=(-4.822564e-02, -1.937633e+00)
Closed curve detected at step 159

npoints curve = 159
>>
The results are plotted in Figure 18 using the standard plot function cpl where the fourth
argument is used to select the fourth and first components of the solution which are the
parameter q2 and the coordinate x.
We not that the BT points are detected twice, namely on both the fold curve and the Hopf
curve. Also, the Hopf curve is a closed curve but only the part between the two BT points that
intersects the equilibrium curve consists of Hopf points. The other part consists of neutral
saddle equilibria which have two real eigenvalues with sum zero.
77
0.16
0.14
0.12 BT
BT
0.1
0.08 H
x
GH
0.06 LP
0.04 CP
LP
0.02 GH H BT
BT
0
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2
q2
Figure 18: Computed equilibrium, fold and Hopf curves in the catalytic oscillator model (60)
8.3 Period Doubling

The Period Doubling bifurcation curve is defined by the following system
 du

 dτ − T f (u, α) =0
 u(0) − u(1) =0
R1 (64)

 hu(t), u̇old (t)idt =0
 0
G(u, T, α) =0
this is exactly the system defining limit cycles but with one extra constraint G(u, T, α) = 0
where G(u, T, α) is defined as the solution component G of the system

 v̇(τ ) − T fu (u, α)v(τ ) + Gϕ(τ ) = 0
v(0) + v(1) = 0 (65)
 R1
0 hψ(τ ), v(τ )idτ = 1
which is exactly the same system as was used to detect the Period Doubling bifurcation.
Instead of using this functional G both systems can be combined in one larger system

 du
 − T f (u, α) =0
 dτ


 u(0) − u(1) =0

 R1
0 hu(t), u̇old (t)idt =0
(66)
 v̇(τ ) − T fu (u, α)v(τ )
 =0



 v(0) + v(1) =0
 R 1 hv (τ ), v(τ )idτ − 1 = 0

0 old
The first method (using system (64) and (65)) is implemented in the curve definition file
perioddoubling. The other method is implemented in perioddoubling2. We will not use
it further and it should not be confused with the continuation of the period doubled curve
that branches off in a PD point on a limit cycle curve.
78
8.3.2 Output of a continuation of period doubling bifurcation points
The x− output is in this case a matrix and each column corresponds to a computed period
doubling limit cycle. The last three components of each column contain the value of the
period T and the values of the two active parameters, in that order.
8.3.3 Bifurcations along a flip curve

In MatCont / CL MatCont there are four generic codim 2 bifurcations that can be
detected along a flip curve:
• Strong 1:2 resonance. We will denote this bifurcation with R2
• Fold - flip, Limit Point - Period Doublingor We will denote this bifurcation with LPPD
• Flip - Neimark-Sacker, denoted as PDNS
• Generalized period doubling point, denoted as GPD
To detect these singularities, we define 4 test functions:
• φ1 = (v1∗ )TW1 LC×M v1M
• φ2 = (ψ1∗ )TW1 (F (u0,1 (t)))C

• φ3 = det ((M ⊙ M ) − In )
• φ4 = c
where c is the coefficient defined in (52), M is the monodromy matrix and ⊙ is the bialternate
product.
The v1 ’s and ψ1 ’s are obtained as follows. For a given ζ ∈ C 1 ([0, 1], Rn ) we consider three
different discretizations:
• ζM ∈ R(N m+1)n the vector of values in the mesh points,
• ζC ∈ RN mn the vector of values in the collocation points,

ζW 1
• ζW = ∈ RN mn × Rn where ζW1 is the vector of values in the collocation points
ζW 2
multiplied with the Gauss - Legendre weights and the lengths of the mesh intervals, and
ζW2 = ζ(0).
Formally we further introduce LC×M which is a structured sparse matrix that converts a
vector ζM of values in the mesh points into a vector ζC of values in the collocation points by
ζC = LC×M ζM . We compute v1M by solving

D − T A(t)
v1M = 0. (67)
δ0 + δ1 disc
P Pm
The normalization of v1M is done by requiring N i=1 j=0 σj h(v1M )(i−1)m+j , (v1M )(i−1)m+j i∆i =
1 where σj is the Gauss-Lagrange quadrature coefficient and ∆i is the length of the i-th in-
terval. By discretization we obtain

∗ T D − T A(t)
(v1W ) = 0.
δ0 + δ1 disc
79
P N Pm R1 ∗
To normalize (v1∗ )W1 we require i=1
∗
j=1 |((v1 )W1 )ij |1 = 1. Then 0 hv1 (t), v1 (t)idt is
approximated by (v1∗ )TW1 LC×M v1M and if this quantity is nonzero, v1W ∗ is rescaled so that
R1 ∗ 1 ∗
0 hv1 (t), v1 (t)idt = 2 . We compute ψ1W by solving

∗ T D − T A(t)
(ψ1W ) =0
δ0 − δ1 disc
∗
P N Pm ∗)
R1 ∗
and normalize ψ1W 1
by requiring i=1 j=1 |((ψ 1 W 1 ) ij | 1 = 1. Then 0 hψ1 (t), F (u0,1 (t))idt
is approximated by (ψ1∗ )TW1 (F (u0,1 (t)))C and if this quantity is nonzero, ψ1W ∗ is rescaled so
R1 ∗ ∗ T
that 0 hψ1 (t), F (u0,1 (t))idt = 1. a1 can be computed as (ψW1 ) (B(t, v1M , v1M ))C . The
computation of (h2,1 )M is done by solving

 (D − T A(t))C×M (h2,1 )M = (B(t; v1M , v1M ))C + 2a1 (F (u0,1 (t)))C , t ∈ [0, 1]
(δ(0) − δ(1))(h2,1 )M = 0,
 ∗ )T L
(ψW 1 C×M (h2,1 )M = 0.
The expression for the normal form coefficient c becomes

c = 13 ((v1W
∗
1
)T (C(t; v1M , T1 v1M , v1M )C + 3(B(t; v1M , (h2,1 )M )C
−6 aT1 (v1W
∗
1
)T (A(t)v1 (t))C ).
The singularity matrix is:  
0 − −
S =  − 0 − . (68)
1 1 0
8.3.4 Period doubling initialization

The natural way to start a Period Doubling bifurcation curve continuation is to start it from a
Period Doubling bifurcation point on a limit cycle curve. This can be done using the following
command: [x0,v0]=init PD PD(@odefile, x, s, ap, ntst, ncol). x should be the x as
returned by the previous limit cycle continuation. s is the special point structure of the
detected Period Doubling point on the limit cycle curve. odefile specifies the ode-file to
be used. ap is the active parameter and ntst and ncol are again the number of mesh and
collocation points for the discretization.
MatCont provides four other initializers to start the continuation of PD curves from a
codim2 bifurcation of limit cycles, namely init GPD PD.m, init LPPD PD.m, init PDNS PD.m
and init R2 PD.m. These initializers are added for ease of use; they only refer back to
init PD PD.m.
We note that MatCont also provides an initializer init PD PD2.m which is merely an
alternative to init PD PD.m; it is not to be confused with the initializer init PD LC.m which
initializes the computation of the period doubled curve from a PD point.
8.3.5 Example
For this example the adaptive control system

 ẋ = y
ẏ = z (69)

ż = −αz − βy − x + x2
80
from §7.7 is used again. The starting vector x0 is calculated from the Period Doubling
bifurcation detected in the limit cycle example (section 7.7) using init PD PD. Continuation
is done using a call to the standard continuer with perioddoubling as curve definition file.
The results are plotted using the standard plot function cpl where the fourth argument is
used to select the 245th and 246th component of the solution which are the parameters α and
β. The results can be seen in Figure 19. The labels of the plot are added manually. It can
be tested by the command testadapt3.
>> [x0,v0]=init_EP_EP(@adaptx,[0;0;0],[-10;1],[1]);
first point found
label = H , x = ( 0.000000 0.000000 0.000000 1.000002 )

npoints curve = 300
>> x1=x(1:3,s(2).index);p=[x(end,s(2).index);1];
>> [x0,v0]=init_H_LC(@adaptx,x1,p,[1],1e-6,20,4);
>> opt = contset(opt,’MaxNumPoints’,200);
>> opt = contset(opt,’Multipliers’,1);
>> opt = contset(opt,’Adapt’,1);
first point found
Neutral Saddle Cycle (period = 6.433818e+00, parameter = 1.895459e-08)

npoints curve = 200
>> plotcycle(xlc,vlc,slc,[size(xlc,1) 1 2]);
>> [x0,v0]=init_PD_PD(@adaptx,xlc,slc(4),[1 2],20,4);
>> [xpd,vpd,spd,hpd,fpd]=cont(@perioddoubling,x0,v0,opt);
first point found
Resonance 1:2 (period = 4.841835e+00, parameters = 5.317604e-09, 1.698711e+00)
(a,b)=(-7.330657e-02, 5.220090e-09)
Resonance 1:2 (period = 9.058318e+00, parameters = -3.045539e-07, 6.782783e-01)
(a,b)=(-1.060571e+01, -3.537444e-04)
81
1.8
R2
1.6
1.4
1.2
0.8
R2
0.6
0.4
-1.5 -1 -0.5 0 0.5 1 1.5
Figure 19: Computed Period Doubling curve

npoints curve = 300
>> cpl(xpd,vpd,spd,[245 246]);
We note that xpd is a 246 × 300 matrix; each column corresponds to a computed period
doubling limit cycle and gives the coordinates of all points of the fine mesh, i.e. 243 =
(20 × 4 + 1) × 3 values, plus the period T as the 244 − th component and the values of the
two active parameter α, β as the 245−th and 246−the components.
8.4 Continuation of fold bifurcation of limit cycles

A Fold bifurcation of limit cycles (Limit Point of Cycles, LPC) generically corresponds to a
turning point of a curve of limit cycles. It can be characterized by adding an extra constraint
G = 0 to (47) where G is the Fold test function. The complete BVP defining a LPC point
using the minimal extended system is
 du

 dt − T f (u, α) =0
 u(0) − u(1) =0
R1 (70)

 hu(t), u̇old (t)idt = 0
 0
G[u, T, α] =0
82
where G is defined by requiring
 
  0
v  0 
N1  S  =  
 0 . (71)
G
1
Here v is a function, S and G are scalars and
 
D − T fu (u(t), α) − f (u(t), α) w01
 δ1 − δ0 0 w02 
N1 = 
 (72)
Intf (u(·),α) 0 w03 
Intv01 v02 0
where the bordering functions v01 , w01 , vector w02 and scalars v02 and w03 are chosen so that
N 1 is nonsingular [15]. This method (using system (71) and (72)) is implemented in the curve
definition file limitpointcycle. The discretization is done using orthogonal collocation in
the same way as it was done for limit cycles.
8.4.2 Bifurcations along a fold of cycles curve

In CL MatCont there are five generic codim 2 bifurcations that can be detected along a
fold of cycles curve:
• Branch point. We will denote this bifurcation with BPC
• Strong 1:1 resonance. We will denote this bifurcation with R1
• Cusp point, denoted by CPC
• Fold - flip, Limit Point - Period Doublingor We will denote this bifurcation by LPPD
• Fold - Neimark-Sacker, denoted by LPNS
A Generalized Hopf (GH) marks the end (or start) of an LPC curve.
To detect the generic singularities, we first define bp+5 test functions, where bp is the number
of branch parameters:
• ψi = wT fβi (u, α), i ∈ (1, ..., bp)
• ψbp+1 = (ϕ∗1 )TW1 LC×M v1M
• ψbp+2 = b (normal form coefficient)
• φbp+3 = det (M + In )
• ψbp+4 = det ((M 2 ⊙ M 2) − In )
In the ψi expressions w is the vector computed in (57) and βi (branch parameter) is a
component of α.
In the second expression ψbp+1 , we compute v1M by solving
   
D − T A(t) T F (u0,1 (t))
 δ0 − δ1  v1M =  0  . (73)
R1
0 F (u0,1 (t))LC×M disc
0 disc
83
By discretization we obtain

D − T A(t)
(ϕ∗1W )T = 0.
δ0 − δ1 disc
P Pm R1 ∗
To normalize (ϕ∗1 )W1 we require N i=1
∗
j=1 ((v1 )W1 )(i−1)m+j 1 = 1. Then 0 hϕ1 (t), v1 (t)idt
is approximated by (ϕ∗1 )TW1 LC×M v1M and if this quantity is nonzero, ϕ∗1W is rescaled so that
R1 ∗
0 hϕ1 (t), v1 (t)idt = 1.
So the third expression for the normal form coefficient b becomes
1
b = ((ϕ∗1W1 )T ((B(t; v1M , v1M )C + 2(A(t)v1 (t))C )).
2
In the fourth expression, M is the monodromy matrix.
In the fifth expression, M 2 is the (n − 2) × (n − 2) matrix that restricts the n × n matrix
M to the space orthogonal to the two-dimensional left eigenspace of the two multipliers that
are closest to 1.
The number of branch parameters is not fixed. If the number of branch parameters is 3
then this matrix has three more rows and columns. This singularity matrix is automatically
extended:  
0 − − − −
 − 0 − − − 
 
S=  − − 0 − − .

 − − − 0 − 
− − − 1 0
8.4.3 Fold of cycles initialization

One way to start a Fold bifurcation of cycles continuation supported in the current ver-
sion is to start it from a fold bifurcation point (LPC) on a limit cycle curve. This can be
done using the following command: [x0,v0]=init LPC LPC(@odefile, x, s, ap, ntst,
ncol(, bp)optional ). x should be the x as returned by the previous limit cycle continuation. s
is the special point structure of the detected Fold bifurcation point on the limit cycle curve.
odefile specifies the ode-file to be used. ap is the active parameter and ntst and ncol
are again the number of mesh and collocation points for the discretization. bp are the op-
tional indices of the branch parameters. It also works without entering a value for this field:
[x0,v0]=init LPC LPC(@odefile, x, s, ap, ntst, ncol).
MatCont provides four other initializers which allow to start a continuation of folds of
limit cycles by starting from an codim2 bifurcation of limit cycles. These are init CPC LPC.m,
init LPNS LPC.m, init LPPD LPC.m and init R1 LPC.m. These initializers are introduced for
ease of use since they refer back to init LPC LPC.m.
A more interesting and indeed nontrivial initializer is init GH LPC.m. Computational
methods to switch to nonhyperbolic cycles from codim 2 bifurcations of equilibria are discussed
in [29].
84
8.4.4 Example: the fast Morris-Lecar equations
We consider the following system

v̇ = y − 0.5(v + 0.5) − 2w(v + 0.7) − m∞ (v − 1)
(74)
ẇ = 1.15(w∞ − w)τ
where m∞ (v) = (1 + tanh((v + 0.01)/0.15))/2, w∞ (v) = (1 + tanh((v − z)/0.145))/2 and

τ = cosh((v − 0.1)/0.29). Here v and w are the state variables and y and z are the parameters.
This is a modification of the fast subsystem of the Morris-Lecar equations studied in [37],[38];
the Morris-Lecar equations were introduced in [33] as a model for the electrical activity in
the barnacle giant muscle fiber. In our model y corresponds to the slow variable in the fast
Morris-Lecar equations; z is the potential for which w∞ (z) = 12 .
We find a stable equilibrium (EP) for y = 0.110472 and z = 0.1 at (0.04722, 0.32564) by
time integration. We continue this equilibrium with free parameter y for decreasing values of
y by running testEquilMLfast.m:
p=[0.11047;0.1];ap1=[1];
[x0,v0]=init_EP_EP(@MLfast,[0.047222;0.32564],p,ap1);
opt=contset;
opt=contset(opt,’MinStepsize’,0.00001);
opt=contset(opt,’MaxStepsize’,0.01);
cpl(x,v,s,[3 1]);
>> testEquilMLfast
first point found
label = H , x = ( 0.036756 0.294770 0.075658 )
label = LP, x = ( -0.033738 0.136501 -0.020727 )
a=-1.036700e+01
Neutral saddle
label = H , x = ( -0.119894 0.045956 0.033207 )
label = LP, x = ( -0.244914 0.008514 0.083257 )
a=-2.697425e+00

npoints curve = 65
We find a Hopf (H) bifurcation point at y = 0.075659, two limit points (LP) at y =
−0.020727 and y = 0.083257 and a neutral saddle (H) at y = 0.033207. There are stable
equilibria before the first H point and after the second LP point and unstable equilibria
85
0.05 H
0
LP
-0.05
-0.1
H
-0.15
-0.2
LP
-0.25
-0.3
-0.35
-0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12
Figure 20: Computed equilibrium curve in the MLfast model
between the first H point and the second LP point. The Lyapunov coefficient in the first
Hopf point l1 = 8.234573 is positive which means that the periodic orbits are born unstable.
The results are plotted using the plot function cpl, cf. §3.9 where the fourth argument is
used to select the third and first component of the solution which are the parameter y and
the coordinate v. The results can be seen in Figure 20.
The detected Hopf point is used to start a limit cycle continuation. We choose N = 30
test intervals and m = 4 collocation points for the discretization.
testEquilMLfast
x1=x(1:2,s(2).index);p=[x(end,s(2).index);0.1];
[x0,v0]=init_H_LC(@MLfast,x1,p,ap1,0.0001,30,4);
opt=contset;
opt=contset(opt,’IgnoreSingularity’,1);
[x2,v2,s2,h2,f2]=cont(@limitcycle,x0,v0,opt);
plotcycle(x2,v2,s2,[1 2]);
>> testLCMLfast
first point found
label = H , x = ( 0.036756 0.294770 0.075658 )
label = LP, x = ( -0.033738 0.136501 -0.020727 )
a=-1.036700e+01
label = H , x = ( -0.119894 0.045956 0.033207 )
label = LP, x = ( -0.244914 0.008514 0.083257 )
86
0.4
0.35
0.3
0.25
w
0.2
0.15
0.1
0.05
−0.1 −0.05 0 0.05 0.1
v
Figure 21: Computed limit cycle curve started from a Hopf bifurcation
a=-2.697425e+00

npoints curve = 65
first point found

npoints curve = 50
The periodic orbit is initially unstable. We detect a limit point of cycles LPC at y =
0.084569. At this point the stability is gained. Afterwards the stability is preserved but the
period tends to infinity and the periodic orbits end in a homoclinic orbit. The results can be
seen in Figure 21.
We now compute a curve of fold bifurcations of limit cycles. The starting vector x0 is cal-
culated from the LPC on the previously computed curve of limit cycles, using init LPC LPC.
Continuation is done using a call to the standard continuer with limitpointcycle as curve
definition file. We free both y and z to continue the LPC curve through this LPC point. The
computations are done by executing the script testLPCMLfast.m:
testEquilMLfast
87
x1=x(1:2,s(2).index);p=[x(end,s(2).index);0.1];
[x0,v0]=init_H_LC(@MLfast,x1,p,ap1,0.0001,30,4);
opt=contset;
opt=contset(opt,’IgnoreSingularity’,1);
opt=contset(opt,’FunTolerance’,0.0000001);
opt=contset(opt,’VarTolerance’,0.0000001);
[x2,v2,s2,h2,f2]=cont(@limitcycle,x0,v0,opt);
[x0,v0]=init_LPC_LPC(@MLfast,x2,s2(2),[1 2],30,4);
opt=contset;
opt=contset(opt,’FunTolerance’,0.0001);
opt=contset(opt,’VarTolerance’,0.0001);
[x3,v3,s3,h3,f3]=cont(@limitpointcycle,x0,v0,opt);
plotcycle(x3,v3,s3,[1 2]);
>> testLPCMLfast
first point found
label = H , x = ( 0.036756 0.294770 0.075658 )
label = LP, x = ( -0.033738 0.136501 -0.020727 )
a=-1.036700e+01
Neutral saddle
label = H , x = ( -0.119894 0.045956 0.033207 )
label = LP, x = ( -0.244914 0.008514 0.083257 )
a=-2.697425e+00

npoints curve = 65
first point found

npoints curve = 50
first point found

npoints curve = 30
88
0.4
0.35
0.3
w
0.25
0.2
0.15
-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12
v
Figure 22: Computed fold of limit cycles curve started from a fold bifurcation of limitcycles
The results are plotted using the standard plot function plotcycle where the fourth
argument is used to select the coordinates. The results can be seen in Figure 22. We note
that it shrinks to a single point. The labels of the plot are added manually .
8.5 Continuation of torus bifurcation of limit cycles

A torus bifurcation of limit cycles (Neimark-Sacker, NS) generically corresponds to a bifur-
cation to an invariant torus, on which the flow contains periodic or quasi-periodic motion.
It can be characterized by adding an extra constraint G = 0 to (47) where G is the torus
test function which has four components from which two are selected. The complete BVP
defining a NS point using a minimally extended system is
 dx

 dt − T f (x, α) =0
 x(0) − x(1) =0
R1 (75)

 hx(t), ẋold (t)idt = 0
 0
G[x, T, α, κ] =0
where
G11 G12
G=
G21 G22
89
is defined by requiring  
  0 0
v1 v2  0 0 
N 3  G11 G12  = 

. (76)
1 0 
G21 G22
0 1
Here v 1 and v 2 are functions and G11 , G12 , G21 and G22 are scalars and
 
D − T fx (x(·), α) w11 w12
 δ0 − 2κδ1 + δ2 w21 w22 
N3 = 
 (77)
Intv01 0 0 
Intv02 0 0
where the bordering functions v01 , v02 , w11 , w12 , vectors w21 and w22 are chosen so that N 3 is
nonsingular [15]. We note that an additional variable κ is introduced in (75). This method
(using system (76) and (77)) is implemented in the curve definition file neimarksacker.m.
The discretization is done using orthogonal collocation over the interval [0 2]. The additional
variable κ is introduced as the last continuation variable of neimarksacker.m, after the state
variables, the period and the two active system parameters.
8.5.2 Bifurcations along a Neimark-Sacker curve

In continuous-time systems there are eight generic codim 2 bifurcations that can be detected
along a torus curve:
• 1:1 resonance. We will denote this bifurcation by R1
• 2:1 resonance point, denoted by R2
• Fold-Neimarksacker point, denoted by LPNS
• Chenciner point, denoted by CH.
• Flip-Neimarksacker point, denoted by PDNS
• Double Neimarksacker bifurcation point, denoted by NSNS

To detect these singularities, we first define 6 test functions:
• ψ1 = κ − 1 (cf. formula (39))
• ψ2 = κ + 1
• ψ3 = κ − 1/2
• ψ4 = κ
• ψ5 = (v1∗ )TW1 LC×M v1M
• ψ6 = Re(d)
90
• ψ7 = det (M + In )
• ψ8 = det ((M 2 ⊙ M 2) − In )
where v1M is computed by solving

D − T A(t) + iθI
v1M = 0. (78)
δ0 − δ1 D
P Pm
The normalization of v1M is done by requiring N i=1 j=0 σj h(v1M )ij , (v1M )ij iti = 1 where
σj is the Gauss-Lagrange quadrature coefficient. By discretization we obtain

∗ H D − T A(t) + iθ
(v1W ) = 0.
δ0 + δ1 disc
P N Pm R1 ∗
To normalize (v1∗ )W1 we require i=1
∗
j=1 |((v1 )W1 )ij |1 = 1. Then 0 hv1 (t), v1 (t)idt is
approximated by (v1∗ )TW1 LC×M v1M and if this quantity is nonzero, v1W ∗ is rescaled so that
R1 ∗ ∗
0 hv1 (t), v1 (t)idt = 1. We compute ϕ1W by solving

∗ T D − T A(t)
(ϕ1W ) =0
δ0 − δ1 disc
P Pm R1 ∗
and normalize ϕ∗1W1 by requiring N i=1
∗
j=1 |((ϕ1 )W1 )ij |1 = 1. Then 0 hϕ1 (t), F (u0,1 (t))idt
is approximated by (ϕ∗1 )TW1 (F (u0,1 (t)))C and if this quantity is nonzero, ϕ∗1W is rescaled so
R1
that 0 hϕ∗1 (t), F (u0,1 (t))idt = 1. We compute h20,1M by solving

D − T A(t) + 2iθI B(t, v1M (t), v1M (t))
h20,1M = .
δ0 − δ1 disc
0
a1 can be computed as (ϕ∗W1 )T (B(t, v1M , v 1M ))C .

The computation of (h11,1 )M is done by solving
   
(D − T A(t))C×M B(t; v1M , v 1M ))C − a1 (F (u0,1 (t)))C
 δ(0) − δ(1)  (h11,1 )M =  0 
∗ T
(ϕW1 ) LC×M 0
The expression for the normal form coefficient d becomes
d = 12 ((v1W
∗
1
)T , (B(t; h11,1M , v1M ))C + (B(t; h20,1M , v 1M ))C + 1
T (C(t; v1M , v1M , v 1M ))C )
− aT1 (v1W
∗
1
)T (A(t)v1 (t))C + iaT12θ .
In the 7th test function, M is the monodromy matrix.
In the 8th test function, M 2 is the (n − 2) × (n − 2) matrix which restricts the n × n matrix
M to the subspace orthogonal to the two-dimensional left eigenspace of the Neimark-Sacker
eigenvalues.
The singularity matrix is:
 
0 − − − − −
 − 0 − − − − 
 
 − − 0 − − − 
S=  . (79)
− − − 0 − − 
 
 − − − − 0 − 
− − − − 1 0
91
8.5.3 Torus bifurcation initialization
One way to start a continuation of torus bifurcations of cycles supported in the current version
is to start it from a torus bifurcation point or neimarksacker point (NS) on a limit cycle curve.
This can be done using the following command: [x0,v0]=init NS NS(@odefile, x, s, ap,
ntst, ncol). x should be the x as returned by the previous limit cycle continuation. s is
the special point structure of the detected torus bifurcation point on the limit cycle curve.
odefile specifies the ode-file to be used. ap is the array containing the two active parameters
and ntst and ncol are again the number of mesh and collocation points for the discretization.
MatCont provides seven other initializers to start the continuation of torus bifurcations
from codim2 bifurcations of limit cycles. These all have the form init XYX NS where XYX is
one of {CH,LP,PD,R1,R2,R3,R4}. They are introduced for ease of use since they all refer
back to init NS NS.m.
More interesting and indeed nontrivial initializers are init HH NS1.m, init HH NS2.m
and init ZH NS. Indeed, two different torus bifurcation curves can generically cross a HH
point. Computational methods to switch to nonhyperbolic cycles from codim 2 bifurcations
of equilibria are discussed in [29].
8.5.4 Example: an autonomous electronic circuit

We consider the following model of an autonomous electronic circuit [18] where x, y and z are
state variables and β, γ, r, a3 , b3 , ν are parameters :

 ẋ = (−(β + ν)x + βy − a3 x3 + b3 (y − x)3 )/r
ẏ = βx − (β + γ)y − z − b3 (y − x)3 (80)

ż = y
A torus bifurcation in this system is described in the manual [14]. It is found by starting
an equilibrium curve from the trival equilibrium point (x = y = z = 0) at β = 0.5, γ = −0.6,
r = 0.6, a3 = 0.32858, b3 = 0.93358, ν = −0.9. The free parameter is ν and the branch is the
trivial one (x = y = z = 0). On this branch a Hopf bifurcation is detected at ν = −0.58934.
On the emerging branch of limit cycles a branch point of limit cycles is found; by continuing
the newly found branch one detects a torus bifurcation of periodic orbits.
We proceed in a somewhat different way; to avoid the branch point of periodic orbits we
start with a slightly perturbed system where the last equation of (80) is replaced by
ż = y + ǫ
where ǫ is a new (artificial) parameter with initially ǫ = 0.001. The perturbed system is
introduced in MatCont in the odefile torBPC.m in the directory Testruns/TestSystems.
It contains a user function ǫ so that the zeroes of the user function correspond to solutions
of the unperturbed system (80).
The trivial solution is replaced by the equilibrium solution x = 0.00125, y = −0.001,
z = 0.00052502 of the perturbed system torBPC.m and we compute a branch of equilibria
with free parameter ν in the now standard way (note that ν is the sixth parameter, ǫ is the
seventh):
>> p=[0.5;-0.6;0.6;0.32858;0.93358;-0.9;0.001];
>> [x0,v0]=init_EP_EP(@torBPC,[0.00125;-0.001;0.00052502],p,[6]);
92
>> opt=contset; opt= contset(opt,’Singularities’,1);
first point found
label = H , x = ( 0.005604 -0.001000 0.002702 -0.589286 )

npoints curve = 10
A Hopf point is found for ν = −0.58928 for the state values x = 0.0056037, y = −0.001,
z = 0.0027021. We start a curve of periodic orbits from this Hopf point, using N = 25 test
intervals and m = 4 collocation points and choosing ν as the free parameter. The results can
be seen in Figure 23.
>> x1=x(1:3,s(2).index);p(6)= x(end,s(2).index);

>> [x0,v0]=init_H_LC(@torBPC,x1,p,[6],0.0001,25,4);
>> opt=contset(opt,’Multipliers’,1);
>> [x,v,s,h,f]=cont(@limitcycle,x0,v0,opt);
first point found

npoints curve = 50
>> plotcycle(x,v,s,[1 2]);
The previous computations are done by running the script testtorBPC1, including drawing
Figure 23 in which the Hopf point and the three bifurcations of limit cycles are clearly visible.
The axis labels were added manually.
In particular, we detect a torus bifurcation point at ν = −0.59575. To recover the torus
bifurcation point of (80) we continue the torus bifurcation in two parameters ν, ǫ. With the
aid of the user function ǫ we will now locate a NS bifurcation of (80). The starting vector x1
is calculated from the NS on this branch using init NS NS. Continuation is done using a call
to the standard continuer with neimarksacker as curve definition file:
p=[0.5;-0.6;0.6;0.32858;0.93358;-0.9;0.001];
x=[0.00125;-0.001;0.00052502];
[x0,v0]=init_EP_EP(@torBPC,x,p,[6]);
opt=contset;
93
0.1
0.05
y
-0.05
-0.1
-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35

x
Figure 23: Computed limit cycle curve started from a Hopf bifurcation
x1=x(1:3,s(2).index);
p(6)=x(end,s(2).index);
[x0,v0]=init_H_LC(@torBPC,x1,p,[6],0.0001,25,4);
opt=contset;
opt=contset(opt,’Multipliers’,1);
[x1,v1]=init_NS_NS(@torBPC,x,s(3),[6 7],25,4);
opt=contset;
opt=contset(opt,’VarTolerance’,1e-4);
opt=contset(opt,’FunTolerance’,1e-4);
opt=contset(opt,’Userfunctions’,1);
UserInfo.name=’epsilon0’;
UserInfo.state=1;
UserInfo.label=’E0’;
opt=contset(opt,’UserfunctionsInfo’,UserInfo);
[xns1,vns1,sns1,hns1,fns1]=cont(@neimarksacker,x1,v1,opt);
plotcycle(xns1,vns1,sns1,[1 2]);
This continuation is done by running the script testtorBPC2. The output is the following:
94
>> testtorBPC2
first point found
label = H , x = ( 0.005604 -0.001000 0.002702 -0.589286 )

npoints curve = 10
first point found

npoints curve = 50
first point found
label = E0, x = ( 0.046835 0.141698 0.046209 ...
... 0.141698 0.046209 8.793572 -0.591612 -0.000006 0.822406 )

npoints curve = 16
We note that the point with label E0 is a torus bifurcation point of (80) with ν = −0.591612,
ǫ = −0.000006 ≈ 0, and period 0.822406. The orbits are shown in Figure 24. The axis labels
were added manually.
Finally, we continue numerically the NS orbit with two free parameters β, ν and find,
interestingly, that the NS orbit shrinks to a single point. The results are plotted using the
standard plot function plotcycle where the fourth argument is used to select the coordinates.
The plot is shown in Figure 25 (the axis labels were added manually). This computation is
executed by running testtorBPC3:
p=[0.5;-0.6;0.6;0.32858;0.93358;-0.9;0.001];
x=[0.00125;-0.001;0.00052502];
[x0,v0]=init_EP_EP(@torBPC,x,p,[6]);
opt=contset;
x1=x(1:3,s(2).index);
p(6)=x(end,s(2).index);
[x0,v0]=init_H_LC(@torBPC,x1,p,[6],0.0001,25,4);
opt=contset;
95
0.1
0.05
y
-0.05
-0.1
-0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2

x
Figure 24: Computed torus bifurcation curve started from a torus bifurcation of limit cycles
opt=contset(opt,’Multipliers’,1);
[x2,v2]=init_NS_NS(@torBPC,x,s(3),[1 6],50,4);
opt=contset; opt=contset(opt,’VarTolerance’,1e-4);
[xns2,vns2,sns2,hns,fns]=cont(@neimarksacker,x2,v2,opt);
plotcycle(xns2,vns2,sns2,[1 2]);
>> testtorBPC3
first point found
label = H , x = ( 0.005604 -0.001000 0.002702 -0.589286 )

npoints curve = 10
first point found
96
0.1
0.05
y
-0.05
-0.1
-0.05 0 0.05 0.1 0.15 0.2

x
Figure 25: Computed torus bifurcation curve that shrinks to a single point


npoints curve = 50
first point found

npoints curve = 50
Here xns2 is a 607 × 50 matrix. The last four components of each column are, in that order,
the period of the orbit, the values of the two free parameters, and the value of the additional
variable κ that is introduced in (75).
97
9 Continuation of codim 2 bifurcations
9.1 Branch Point Continuation
Although the predecessor CONTENT [26] of MatCont included the continuation of codi-
mension 2 bifurcations of equilibria in three free parameters, this is so far not included in
MatCont. Though the initialization and continuation itself would not be too difficult, the
computation and interpretation of the normal form coefficients of the codim3 bifurcations
would be unfeasible. Nevertheless, MatCont includes the continuation of branch points of
equilibria and limit cycles in three parameters since this appears to be useful in some studies
of models with equivariance.

In the toolbox branch point curves are computed by minimally extended defining systems, cf.
[22], §4.1.2. The branch point curve is defined by the following system

 f (u, α) = 0,
g1 (u, α) = 0, (81)

g2 (u, α) = 0,
where (u, α) ∈ Rn+2 , while g1 and g2 are obtained by solving
   
v11 v21 0n 0n
N 4  v12 v22  =  1 0  . (82)
g1 g2 0 1
Here v11 and v21 are functions and v12 , v22 , g1 and g2 are scalars and
 
fu (u, α) fβ (u, α) w01
N 4 =  v011T v012T 0 
v021T v022T 0
where the bordering functions v011 , v021 , w01 and scalars v012 , v022 are chosen so that N 4 is non-
singular. This method is implemented in the curve definition file branchpoint.m.
9.1.2 Bifurcations
In the current version no bifurcations are detected.
9.1.3 Branch Point initialization

The only way to start a continuation of branch points supported in the current version is to
start it from a Branch point detected on an equilibrium curve or on a fold curve. This can
be done using the following statement: [x0,v0]=init BP BP(@odefile, x, p, ap, bp).
This routine stores its information in a global stucture bpds. The result of init BP BP
contains a vector x0 with the state variables and the three active parameters and a vector
v0 that is empty. Here odefile is the ode-file to be used, x is a vector of state variables
containing the values of the state variables returned by a previous equilibrium curve or fold
curve continuation, p is the vector containing the current values of the parameters, ap is the
vector containing the indices of the 3 active parameters and bp is the index of the branch
parameter.
98
9.1.4 Example
For this example we use a model of a continuous stirred tank reactor (cstr.m):

ẋ = α3 − (1 + λ)x + λα1 /(1 + λα2 ∗ e−α4 x/(1+x) ) (83)
where α1 = 10 − 9 ∗ β + γ, α2 = 10 − 9 ∗ β and α3 = −0.9 + 0.4 ∗ β. The model is coded in

such a way that the parameters are λ, β, γ, α4 , in that order.
It is easily seen that x = −0.9 is an equilibrium point of the system for the choice
(0, 0, 0, 3) of the parameters. From this we can start an equilbrium continuation with λ (the
first parameter) free.
p=[0;0;0;3];ap1=[1];
[x0,v0]=init_EP_EP(@cstr,[-0.9],p,ap1);
opt=contset;
opt=contset(opt,’VarTolerance’,1e-3);
first point found
label = LP, x = ( -0.143564 1.250669 )
a=1.550147e+000
label = LP, x = ( 0.393180 0.377651 )
a=-7.370472e-001
cpl(x,v,s,[2,1]);
The results are plotted using the plot function cpl where the fourth argument is used to
select the second and first components of the solution which are the parameter λ and the
coordinate x. The resulting curve is a part of Figure 26. These computations can be done by
running the script cstr1.
We start a fold continuation from the second LP detected on the previous equilibrium
curve; λ and β are free in this run.
x1=x(1,s(3).index);
[x0,v0]=init_LP_LP(@cstr,x1,p,[1 2],[1 2 3 4]);
first point found
label = BP1, x = ( 2.018621 0.581081 -4.709219 )
label = CP , x = ( 0.259553 1.968966 -0.090655 )
c=-8.847089e-001
label = BP1, x = ( 0.030643 1.772454 -0.127542 )
label = BP4, x = ( -0.000009 1.707401 -0.124964 )
label = CP , x = ( -0.173872 0.405524 0.608093 )
99
2.5
BP1
2
1.5
0.5 LP
CP
BP4 BP1
BP4
0 CP LP
−0.5
−1
−0.5 0 0.5 1 1.5 2
Figure 26: Computed fold curve
c=-2.263137e+000
label = BP4, x = ( -0.000000 0.421692 0.528995 )

npoints curve = 164
hold on;
cpl(x2,v2,s2,[2,1]);
These computations can be done by running the script cstr2. The results are plotted
using the standard plot function cpl where the fourth argument is used to select the second
and first components of the solution which are the parameter λ and the coordinate x. The
results can be seen in Figure 26.
Finally, we continue numerically the BP curves with three free parameters λ, β and γ.
The BP curves are started respectively from the first BP1 point (λ is the branch parameter)
and the first BP4 point (α4 is the branch parameter) detected on the previous fold curve.
The results are plotted using the standard plot function cpl where the fourth argument is
used to select the coordinates. A graphical representation of this phenomenon is shown in
Figure 27. In the latter λ is plotted versus x. The labels of the plot are added manually .
x1=x2(1,s2(2).index);
p([1 2])=x2(end-1:end,s2(2).index);
[x0,v0]=init_BP_BP(@cstr,x1,p,[1 2 3],1);
[x3,v3,s3,h3,f3]=cont(@branchpoint,x0,[],opt);
first point found

npoints curve = 300
100
2.5
BP1
2
1.5
0.5 LP
CP
BP1
0 BP4 LP BP4
CP
−0.5
−1
−0.5 0 0.5 1 1.5 2 2.5 3 3.5
Figure 27: Computed BP curves started from Branch Points detected on a fold curve.
hold on;
cpl(x3,v3,s3,[2,1]);
x1=x2(1,s2(5).index);
p([1 2])=x2(end-1:end,s2(5).index);
[x0,v0]=init_BP_BP(@cstr,x1,p,[1 2 3],4);
[x3,v3,s3,h3,f3]=cont(@branchpoint,x0,[],opt);
first point found

npoints curve = 300
hold on;
cpl(x3,v3,s3,[2,1]);
These computations can be done by running the script cstr3.
9.2 Branch Point of Cycles Continuation

9.2.1 Mathematical Definition
A BPC can be characterized by adding two extra constraints G1 = 0 and G2 = 0 to (47)
where G1 and G2 are the Branch Point test functions. The complete BVP defining a BPC
101
point using the minimal extended system is
 dx

 dt − T f (x, α) =0
 x(0) − x(1) =0
R1 (84)

 hx(t), ẋold (t)idt =0
 0
G[x, T, α] =0
where
G1
G=
G2
is defined by requiring  
0 0
 0 0 
v1 v2  
N =
 0 0 .
 (85)
G1 G2  
1 0
0 1
Here v1 and v2 are functions, G1 and G2 are scalars and
 
D − T fx (x(·), α) − f (x(·), α) − T fβ (x(·), α) w01
 δ 1 − δ 0 0 0 w02 
 
N =  Intẋold (·) 0 0 w03  (86)

 v11 v12 v13 0 
v21 v22 v23 0
where the bordering operators v11 , v21 , function w01 , vector w02 and scalars v12 , v22 , v13 , v23
and w03 are chosen so that N is nonsingular [15][16].
9.2.2 Bifurcations
In the current version no bifurcations are detected.
9.2.3 Branch Point of Cycles initialization

The only way to start a continuation of branch points supported in the current version is to
start it from a BPC detected on an limitcycle curve or on an LPC curve. This can be done
using the following statement: [x0,v0]=init BPC BPC(@odefile, x, s, ap, ntst, ncol,
bp). This routine stores its information in a global stucture lds. The result of init BP BP
contains a vector x0 with the state variables and the three active parameters and a vector
v0 that is empty. Here odefile is the ode-file to be used, x is a vector of state variables
containing the values of the state variables returned by a previous limit cycle curve or fold
of cycles curve continuation, s is a structure containing the BPC point values returned by a
previous limit cycle curve or fold of cycles curve continuation, ap is the vector containing the
indices of the 3 active parameters and bp is the index of the branch parameter. ntst and
ncol are again the number of mesh and collocation points for the discretization.
102
9.2.4 Example
In this section we discuss a non-generic situation, i.e. a case with a symmetry and a contin-
uation of BPC points that involves two effective parameters and one artificial parameter.
For this example the following model is used:

 ẋ = (−(β + ν)x + βy − a3 x3 + b3 (y − x)3 )/r
ẏ = βx − (β + γ)y − z − b3 (y − x)3 (87)

ż = y
which is the same system as in the torus of cycles continuation. It has a trivial solution branch
x = y = z = 0 for all parameter values. Moreover, it has the Z2 - symmetry x → −x,y → −y.
To compute the branch of BPC points with respect to ν through the BPC point that we will
detect on a limitcycle continued with free parameters ν, β, we need to introduce an additional
free parameter that breaks the symmetry. There are many choices for this; we choose to
introduce a parameter ǫ and extend the system (87) by simply adding a term +ǫ to the first
right-hand-side. For ǫ = 0 this reduces to (87) while for ǫ 6= 0 the symmetry is broken.
We start by computing the trivial branch with fixed parameters γ = −0.6, r = 0.6,
a3 = 0.328578, b3 = 0.933578, β = 0.5,ǫ = 0 and free parameter ν with initially ν = −0.9.
On this branch a Hopf point is detected for ν = −0.58933644 and a branch point of equilibria
for ν = −0.5.
>> p=[0.5;-0.6;0.6;0.32858;0.93358;-0.9;0];
>> [x0,v0]=init_EP_EP(@torBPC,[0;0;0],p,[6]);
first point found
label = H , x = ( 0.000000 0.000000 0.000000 -0.589336 )
label = BP, x = ( 0.000000 0.000000 0.000000 -0.500000 )

npoints curve = 50
These computations can be done by running the script testtorBPC4.

From the Hopf point we start the computation of a curve of limit cycles, using 25 test
intervals and 4 collocation points. This is clearly a branch of symmetric solutions of (87); we
detect one LPC and two BPC, see Fig. 28.
>> x1=x(1:3,s(2).index);p(6)= x(end,s(2).index);ap = 6;

>> [x0,v0]=init_H_LC(@torBPC,x1,p,ap,0.0001,25,4);
>> opt=contset(opt,’Multipliers’,1);
>> opt=contset(opt,’Adapt’,5);
first point found
103
0.2
0.15
0.1
0.05
0
y
-0.05
-0.1
-0.15
-0.2
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
x
Figure 28: Curve of limit cycles with LPC and branch points in the circuit example.

Branch Point cycle(period = 8.689669e+000, parameter = -5.870290e-001)
Neutral saddle

npoints curve = 150
>> plotcycle(xlc,vlc,slc,[1 2]);
We continue the secondary cycle branch passing through the BPC point. From Fig. 29 it
is clear that in the secondary cycle the symmetry is broken.
>> [x1,v1]=init_BPC_LC(@torBPC,xlc,vlc,slc(3),25,4,1e-6);
>> opt=contset(opt,’Backward’,1);
>> [xlc1,vlc1,slc1,hlc1,flc1]=cont(@limitcycle,x1,v1,opt);
first point found
104
0.15
0.1
0.05
0
y
−0.05
−0.1
−0.15
−0.1 0 0.1 0.2 0.3 0.4

x
Figure 29: Asymmetric curve of limit cycles in the circuit example.

npoints curve = 50
>> plotcycle(xlc1,vlc1,slc1,[1 2]);
Using the code for the continuation of generic BPC points with three free parameters
ν, β, ǫ we continue the curve of non-generic BPC points, where ǫ remains close to zero. The
picture in Fig. 30 clearly shows that the symmetry is preserved (the axis labels were added
manually).
>> [x1,v1]=init_BPC_BPC(@torBPC,xlc,slc(3),[1 6 7],25,4,ap);
>> [xbpc,vbpc,sbpc,hbpc,fbpc]=cont(@branchpointcycle,x1,v1,opt);
first point found

npoints curve = 200
>> plotcycle(xbpc,vbpc,sbpc,[1 2]);
105
0.8
0.6
0.4
0.2
0
y
-0.2
-0.4
-0.6
-0.8
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

x
Figure 30: Curve of BPC points in the circuit example.
106
10 Continuation of homoclinic and heteroclinic orbits
10.1 Homoclinic orbits: Mathematical definition
In dynamical systems theory, an orbit corresponding to a solution u(t) is called homoclinic
to the equilibrium point u0 of the dynamical system if u(t) → u0 as t → ±∞. There are
two types of homoclinic orbits with codimension 1, namely homoclinic-to-hyperbolic-saddle
(HHS), whereby u0 is a saddle, saddle-focus or bi-focus, and homoclinic-to-saddle-node (HSN),
whereby u0 is a saddle-node (i.e., exhibits a limit point bifurcation). We recall that auto
has a toolbox for homoclinic continuation, named HomCont [4], [5].
During the continuation, it is necessary to keep track of several eigenspaces of the equi-
librium in each step. To do this in an efficient way, MatCont incorporates the continuation
of invariant subspaces [7] into the defining system.
10.1.1 Homoclinic-to-Hyperbolic-Saddle Orbits

To continue HHS orbits in two free parameters, we use an extended defining system that
consists of several parts.
First, the infinite time interval is truncated, so that instead of [−∞, +∞] we use [−T, +T ],
which is scaled to [0, 1] and divided into mesh-intervals. The mesh is nonuniform and adap-
tive. Each mesh interval is further subdivided by equidistant fine mesh points. Also, each
mesh interval contains a number of collocation points. (This discretization is the same as
that in auto for boundary value problems.) The equation
ẋ(t) − 2T f (x(t), α) = 0, (88)
must be satisfied in each collocation point.
The second part is the equilibrium condition
f (x0 , α) = 0. (89)
Third, there is a so-called phase condition needed for the homoclinic solution, similar to
periodic solutions Z 1
∗
x
ė (t)[x(t) − x
e(t)]dt = 0. (90)
0
Here xe(t) is some initial guess for the solution, typically obtained from the previous continua-
tion step. We note that in the literature another phase condition is also used, see, for example
[13]. However, in the present implementation we employ the condition (90).
Fourth, there are the homoclinic-specific constraints to the solution. For these we need
access to the stable and unstable eigenspaces of the system in the equilibrium point after
each step. It is not efficient to recompute the spaces from scratch in each continuation-step.
Instead, we use the algorithm for continuing invariant subspaces, as described in [7]. This
method adds two small-sized vectors (YS and YU ) to the system variables, from which the
necessary eigenspaces (stable and unstable, respectively) can easily be computed in each step.
107
If Q(0) is an orthogonal matrix whose first m columns form a basis for the invariant sub-
space under consideration in the previous step, and A = fx (x0 , α) is the Jacobian at the new
equilibrium point, then we first compute the so-called Ricatti-blocks, Tij , by the formula

T11 T12
= Q(0)T A Q(0). (91)
T21 T22
If n is the number of state variables, then T11 is of size m × m and T22 is (n − m) × (n − m).
This is done for the stable and unstable eigenspaces separately. Now YS and YU are obtained
from the Ricatti equations
T22U YU − YU T11U + T21U − YU T12U YU = 0,

(92)
T22S YS − YS T11S + T21S − YS T12S YS = 0.
Now we can formulate constraints on the behavior of the solution close to the equilibrium
x0 . The initial vector of the orbit, (x(0) − x0 ), is placed in the unstable eigenspace of the
system in the equilibrium. We express that by the requirement that it is orthogonal to the
orthogonal complement of the unstable eigenspace. Using YU , we can compute the orthogonal
complement of the unstable eigenspace. If QU (0) is the orthogonal matrix from the previous
step, related to the unstable invariant subspace, then a basis for the orthogonal complement
in the new step Q⊥ U (s) is
⊥ −YUT
QU (s) = QU (0) .
I
Note that Q⊥
U (s) is not orthogonal. The full orthogonal matrix Q1U needed for the next step,
is computed separately after each step. The equations to be added to the system are (after
analogous preparatory computations for the stable eigenspace)
Q⊥ T
U (s) (x(0) − x0 ) = 0, (93)
⊥ T
QS (s) (x(1) − x0 ) = 0.
Finally, the distances between x(0) (resp., x(1)) and x0 must be small enough, so that
kx(0) − x0 k − ǫ0 = 0,
(94)
kx(1) − x0 k − ǫ1 = 0.
A system consisting of all equations (88), (89), (90), (92), (93) and (94), is overdetermined.
The basic defining system for the continuation of a HHS orbit in two free parameters consists
of (88), (89), (92), (93), and (94) with fixed ǫ0,1 , so that the phase condition (90) is not used.
The variables in this system are stored in one vector. It contains the values of x(t) in the
fine mesh points including x(0) and x(1), the truncation time T , two free system parameters,
the coordinates of the saddle x0 , and the elements of the matrices YS and YU . Alternatively,
the phase condition (90) can be added if T is kept fixed but ǫ0 and ǫ1 are allowed to vary.
It is also possible to fix T and ǫ0 , say, and allow ǫ1 to vary, again with no phase condition.
Other combinations are also possible, in particular, when the homotopy method [7] is used
to compute a starting homoclinic solution. For more details on the implementation of the
homoclinic continuation we refer to [19].
108
10.1.2 Homoclinic-to-Saddle-Node Orbits
For a homoclinic orbit to a saddle-node equilibrium, the extended defining system undergoes
some small changes. Now (x(0) − x0 ) has to be placed in the center-unstable subspace.
Analogously, (x(1) − x0 ) must be in the center-stable subspace. This again is implemented
by requiring that the vector is orthogonal to the orthogonal complement of the corresponding
space. So the equations (93) themselves do not really change; the changes happen in the
computation of the matrices Q. The defining system now has one equation less than in
the HHS case (ns + nu < n, with ns the dimension of the stable, and nu of the unstable
eigenspace); the number of equations is restored however, by adding the constraint that the
equilibrium must be a saddle-node. For this we use the bordering technique, as described in
section 4.2.1 of [22].
10.2 Bifurcations
During HSN continuation, only one bifurcation is tested for, namely the non-central homoclinic-
to-saddle-node orbit or NCH. This orbit forms the transition between HHS and HSN curves.
The strategy used for detection is taken from HomCont [5].
During HHS continuation, all bifurcations detected in HomCont are also detected in our
implementation. For this, mostly test functions from [5] are used.
Suppose that the eigenvalues of fx (x0 , α0 ) can be ordered according to
ℜ(µns ) ≤ ... ≤ ℜ(µ1 ) < 0 < ℜ(λ1 ) ≤ ... ≤ ℜ(λnu ), (95)
where ℜ() stands for ’real part of’, ns is the number of stable, and nu the number of unstable
eigenvalues. The test functions for the bifurcations are
• Neutral saddle, saddle-focus or bi-focus
ψ = ℜ(µ1 ) + ℜ(λ1 )
If both µ1 and λ1 are real, then it is a neutral saddle, if one is real and one consists of a
pair of complex conjugates, the bifurcation is a saddle-focus, and it is a bi-focus when
both eigenvalues consist of a pair of complex conjugates.
• Double real stable leading eigenvalue

(ℜ(µ1 ) − ℜ(µ2 ))2 , ℑ(µ1 ) = 0
ψ=
−(ℑ(µ1 ) − ℑ(µ2 ))2 , ℑ(µ1 ) 6= 0
• Double real unstable leading eigenvalue

(ℜ(λ1 ) − ℜ(λ2 ))2 , ℑ(λ1 ) = 0
ψ=
−(ℑ(λ1 ) − ℑ(λ2 ))2 , ℑ(λ1 ) 6= 0
• Neutrally-divergent saddle-focus (stable)
ψ = ℜ(µ1 ) + ℜ(µ2 ) + ℜ(λ1 )
109
• Neutrally-divergent saddle-focus (unstable)
ψ = ℜ(µ1 ) + ℜ(λ2 ) + ℜ(λ1 )
• Three leading eigenvalues (stable)
ψ = ℜ(µ1 ) − ℜ(µ3 )
• Three leading eigenvalues (unstable)
ψ = ℜ(λ1 ) + ℜ(λ3 )
• Non-central homoclinic-to-saddle-node
ψ = ℜ(µ1 )
• Shil’nikov-Hopf
ψ = ℜ(λ1 )
• Bogdanov-Takens point
ℜ(µ1 )
ψ=
ℜ(λ1 )
For orbit- and inclination-flip bifurcations, we assume the same ordering of the eigenvalues
of fx (x0 , α0 ) = A(x0 , α0 ) as shown in (95), but also that the leading eigenvalues µ1 and λ1
are unique and real:
ℜ(µns ) ≤ ... ≤ ℜ(µ2 ) < µ1 < 0 < λ1 < ℜ(λ2 ) ≤ ... ≤ ℜ(λnu ) .
Then it is possible to choose normalised eigenvectors ps1 and pu1 of AT (x0 , α0 ) and q1s and q1u
of A(x0 , α0 ) depending smoothly on (x0 , α0 ), which satisfy
AT (x0 , α0 ) ps1 = µ1 ps1 AT (x0 , α0 ) pu1 = λ1 pu1

A(x0 , α0 ) q1s = µ1 q1s A(x0 , α0 ) q1u = λ1 q1u .
The test functions for the orbit-flip bifurcations are then:
• Orbit-flip with respect to the stable manifold
ψ = e−µ1 T < ps1 , x(1) − x0 >
• Orbit-flip with respect to the unstable manifold
ψ = eλ1 T < pu1 , x(0) − x0 >
For the inclination-flip bifurcations, in [5] the following test functions are introduced:
• Inclination-flip with respect to the stable manifold
ψ = e−µ1 T < q1s , φ(0) >
110
• Inclination-flip with respect to the unstable manifold
ψ = eλ1 T < q1u , φ(1) >
where φ (φ ∈ C 1 ([0, 1], Rn )) is the solution to the adjoint system, which can be written as


 φ̇(t) + 2 T AT (x(t), α0 ) φ(t) = 0
 (L )T φ(1) = 0
s
(96)

 (L u )T φ(0) = 0
R
 1 eT e
0 φ (t)[φ(t) − φ(t)]dt = 0
where Ls and Lu are matrices whose columns form bases for the stable and unstable eigenspaces
of A(x0 , α0 ), respectively, and the last condition selects one solution out of the family cφ(t)
for c ∈ R. Lu is equivalent to QU from the mathematical definition of the system, and Ls to
QS . In the homoclinic defining system the orthogonal complements of QS and QU are used;
in the adjoint system for the inclination-flip bifurcation, we use the matrices themselves (or
at least, their transposed versions).
10.3 Homoclinic initialization (HHS)

For homoclinics the same problems occur as with limit cycles, a homoclinic continuation
cannot be done by just calling the continuer as
[x,v,s,h,f]=cont(@homoclinic, x0, v0, opt)
The homoclinic curve file has to know
• which (two) system parameters are active,
• the number of mesh and collocation points to use for the discretization.
Also, some initial information about the state variables is necessary. In the simplest cases
this is an initial cycle x0 with a large period and close to a homoclic orbit, or an already
known homoclinic orbit. This information can be supplied using an initializer:
• [x0,v0]=init LC Hom(@odefile, x, s, p, ap, ntst, ncol, extravec, T, eps0,

eps1)
Calculates an initial homoclinic orbit from a limit cycle with large period. Here odefile
is the ode-file to be used. x and s are here the x and s belonging to the limit cycle
with large period, obtained in a previous continuation. p is the vector containing the
current values of the parameters. ap is the active parameter and ntst and ncol are the
number of mesh and collocation points to be used for the discretization. extravec is
a vector of 3 integers, which are either 0 or 1, and which indicate which of T, eps0,
eps1 are to be variable during the continuation. The vector is 1 for those that should
be variable. This can either be 1 or 2 or the three parameters. T, eps0 and eps1 are
values for these parameters.
111
• [x0,v0]=init Hom Hom(@odefile, x, s, p, ap, ntst, ncol, extravec, T, eps0,
eps1)
Calculates an initial homoclinic orbit from a homoclinic orbit obtained during a previous
continuation. All parameters are similar to the initialisations above.
These initializers return an initial homoclinic orbit x0 as well as its tangent vector v0.
It is also possible to start a homoclinic orbit from a Bogdanov-Takens equilibrium point.
The initialization of homoclinic orbits from a BT point has a long history. We refer in
particular to [30], [31] and [1] where further references can be found. This case is relegated
to the MatCont GUI tutorials.
The homotopy method even allows a start from any equilibrium point. For details and
examples we refer to [8]. This case is also relegated to the MatCont GUI tutorials.
10.4 Homoclinic-to-Saddle-Node initialization (HSN)

For homoclinic-to-saddle-node orbits, the user needs to replace Hom with HSN in the above
two initializers. The homotopy method allows a start from any limit point. For details and
examples we refer to [8]. This case is relegated to the MatCont GUI tutorials.
10.5 CL MatCont: the MLFast example

In §8.4.4 we studied a continuation of limit cycles in the MLfast Morris-Lecar model and
noted that the limit cycles were approaching a homoclinic orbit. We will now approach this
homoclinic even closer, and then start its continuation from the large limit cycle. The result
is shown in Figure 31.
>> init;
>> p=[0.11047;0.1];ap1=[1];
>> [x0,v0]=init_EP_EP(@MLfast,[0.047222;0.32564],p,ap1);
>> opt=contset;opt=contset(opt,’Singularities’,1);
>> opt=contset(opt,’MinStepSize’,0.00001);
>> opt=contset(opt,’MaxStepSize’,0.01);
>> opt=contset(opt,’Backward’,1);
first point found
label = H , x = ( 0.036756 0.294770 0.075659 )
label = LP, x = ( -0.033738 0.136501 -0.020727 )
a=-1.036706e+001
label = H , x = ( -0.119894 0.045956 0.033207 )
Neutral saddle
label = LP, x = ( -0.244915 0.008514 0.083257 )
a=2.697414e+000

npoints curve = 65
112
>> x1=x(1:2,s(2).index);p=[x(end,s(2).index);0.1];
>> [x0,v0]=init_H_LC(@MLfast,x1,p,ap1,0.0001,30,4);
>> opt=contset;
>> opt=contset(opt,’MaxStepSize’,1);
>> opt=contset(opt,’IgnoreSingularity’,1);
>> opt=contset(opt,’Singularities’,1);
>> [x2,v2,s2,h2,f2]=cont(@limitcycle,x0,v0,opt);
first point found

npoints curve = 200
>> p(ap1) = x2(end,end);
>> T = x2(end-1,end)/2;
>> [x0,v0]=init_LC_Hom(@MLfast, x2(:,end), s2(:,end), p, [1 2], 40, 4,...
>> [0 1 1], T, 0.01, 0.01);
>> [xh,vh,sh,hh,fh] = cont(@homoclinic,x0,v0,opt);
first point found
npoints curve = 15
>> plotcycle(xh,vh,sh,[1 2]);
The above computations can be done by running the script homoc1. The picture is pre-
sented in Figure 31. Similar tests can be done by using the testrun testmyml.
10.6 Heteroclinic orbits (Het)

An orbit corresponding to a solution u(t) is called heteroclinic if there exist equilibria u0 , u1 , u0 6=
u1 so that u(t) → u0 as t → −∞ and u(t) → u1 as t → ∞.
Details on the continuation of heteroclinic orbits can be found in [19] and [8]. The routines
for dealing with heteroclinic orbits are in the directory Heteroclinic.
The directories HomotopyHet, HomotopySaddle and HomotopySaddleNode are initializa-
tion directories in which homotopy methods are provided to initialize heteroclinic orbits,
orbits homoclinic to saddle and orbits homoclinicto saddle-node, respectively. Eamples of
their use are relegated to the MatCont GUI tutorials.
113
0.4
0.35
0.3
0.25
w
0.2
0.15
0.1
0.05
−0.15 −0.1 −0.05 0 0.05 0.1

v
Figure 31: Computed curve of homoclinic-to-saddle orbits started from a limit cycle with
large period.
114
A Continuer example: a curve object
In this section a simple example is presented to illustrate the basic use of the continuer. This
example generates a curve g in the (x, y)-plane such that x2 + y 2 = 1. So if the user specifies
a point reasonably close to this curve we get the unit circle. The defining function is
F (x, y) = x2 + y 2 − 1 (97)
In the following listing this curve is implemented. Do note that while there are no options
needed, the curve file must return an option structure (see section 3.3).
curve.m
1 function out = curve
2 %
3 % Curve file of circle
4 %
5
6 out{1} = @curve_func;
7 out{2} = @defaultprocessor;
8 out{3} = @options;
9 out{4} = []; %@jacobian;
10 out{5} = []; %@hessians;
11 out{6} = []; %@testf;
12 out{7} = []; %@userf;
13 out{8} = []; %@process;
14 out{9} = []; %@singmat;
15 out{10} = []; %@locate;
16 out{11} = []; %@init;
17 out{12} = []; %@done;
18 out{13} = @adapt;
19 function f = curve_func(arg)
20 x = arg;
21 f = x(1)^2+x(2)^2-1;
22
23 function varargout = defaultprocessor(varargin)
24 if nargin > 2
25 s = varargin{3};
26 varargout{3} = s;
27 end
28 % no special data
29 varargout{2} = [];
30 % all done succesfully
31 varargout{1} = 0;
32
33 function option = options
34 option = contset;
35
36 function [res,x,v] = adapt(x,v)
37 res=[];
38
curve.m
The file curve.m is stored in the directory Testruns/TestSystems. Starting computations
115
0.8
0.6
0.4
x(2) 0.2
−0.2
−0.4
−0.6
−0.8
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

x(1)
Figure 32: Computed curve of curve.m
at (x, y) = (1, 0), the output in matlab looks like:
>> init;
>> [x,v,s]=cont(@curve,[1;0]);
first point found

npoints curve = 70
The generated curve is plotted in Figure 32 with the command:
>> cpl(x,v,s)
In this case x has dimension 2, so a 2D-plot is drawn with the first component of x (value of
the state variable x) on the x-axis and the second component of x (value of the state variable
y) on the y-axis. The above commands are also executed by running testdrawcurve; the file
testdrawcurve.m is in the directory Testruns.
116
B The Brusselator example: Continuation of a solution to a
boundary value problem in a free parameter
Discretized solutions of PDEs can also be continued in CL MatCont. We illustrate this
by continuing the equilibrium solution to a one-dimensional PDE. The curve type is called
’pde 1’. It is defined in the file pdf 1.m in the directory Equilibrium.
The Brusselator is a system of equations intended to model the Belusov-Zhabotinsky
reaction. This is a system of reaction-diffusion equations that is known to exhibit oscillatory
behavior. The unknowns are the concentrations X(x, t), Y (x, t), A(x, t) and B(x, t) of four
reactants. Here t denotes time and x is a one-dimensional space variable normalized so that
x ∈ [0, 1]. The length L of the reactor is a parameter of the problem. In our simplified setting
A and B are constants.
The system is described by two partial differential equations:
∂X Dx ∂2X
∂t = L2 ∂x2
+ A − (B − 1)X + X 2 Y
∂Y Dy ∂2Y (98)
∂t = L2 ∂x2
+ BX − X 2 Y
with x ∈ [0, 1], t ≥ 0. Here Dx , Dy are the diffusion coefficients of X and Y . At the boundaries
x = 0 and x = 1 Dirichlet conditions are imposed:

X(0, t) = X(1, t) = A
(99)
Y (0, t) = Y (1, t) = B
A
We are interested in equilibrium solutions X(x) and Y (x) to the system and their dependence
on the parameter L.
The approximate equilibrium solution is:

X(x) = A + 2 sin(πx)
(100)
Y (x) = B 1
A − 2 sin(πx)
The initial values of the parameters are: A = 2, B = 4.6, Dx = 0.0016, Dy = 0.08

and L = 0.06. The initial solution (100) is not an equilibrium, but the continuer will try
to converge to an equilibrium close to the initial solution. We use equidistant meshes. To
avoid spurious solutions (solutions that are induced by the discretization but do not actually
correspond to solutions of the undiscretized problem) one can vary the number of mesh points
by setting the parameter N . If the same solution is found for several discretizations, then we
can assume that they correspond to solutions of the continuous problem.
The second order space derivative is approximated using the well-known three-points
2
difference formula: ∂∂xf2 = h12 (fi−1 − 2fi + fi+1 ), where h = N 1+1 , where N is the number of
grid points on which we discretize X and Y . So N is a parameter of the problem and 2N is
the number of state variables (which is not fixed in this case).
The Jacobian is a sparse 5-band matrix. In the ode-file describing the problem the Ja-
cobian is introduced as a sparse matrix. The Hessian is never computed as such but second
order derivatives are computed by finite differences whenever needed. We note that matlab
6.5 or 7 does not provide sparse structures for 3-dimensional arrays.
bruss.m
117
1 function out = bruss
2 %
3 % Odefile of 1-d Brusselator model
4 %
5
6 out{1} = @init;
7 out{2} = @fun_eval;
8 out{3} = @jacobian;
9 out{4} = @jacobianp;
10 out{5} = [];%@hessians;
11 out{6} = [];%@hessiansp;
12 out{7} = [];
13 out{8} = [];
14 out{9} = [];
15
16
17
18 % --------------------------------------------------------------------------
19
20
21 function dfdt = fun_eval(t,y,N,L)
22
23 x = y(1:N);
24 y = y(N+1:2*N);
25
26 A = 2;
27 B = 4.6;
28 Dx = 0.0016;
29 Dy = 0.008;
30 x0 = A; x1 = A;
31 y0 = B/A; y1 = B/A;
32 L2 = L^2;
33 h = 1/(N+1);
34 cx = (Dx/L2)/(h*h);
35 cy = (Dy/L2)/(h*h);
36
37 dxdt = zeros(N,1);
38 dydt = zeros(N,1);
39
40 dxdt(1) = (x0-2*x(1)+x(2))*cx + A - (B+1)*x(1) + x(1)*x(1)*y(1);
41 dxdt(N) = (x(N-1)-2*x(N)+x1)*cx + A - (B+1)*x(N) + x(N)*x(N)*y(N);
42
43 dydt(1) = (y0-2*y(1)+y(2))*cy + B*x(1) - x(1)*x(1)*y(1);
44 dydt(N) = (y(N-1)-2*y(N)+y1)*cy + B*x(N) - x(N)*x(N)*y(N);
45
46 for i=2:N-1
47 dxdt(i) = (x(i-1)-2*x(i)+x(i+1))*cx + A - (B+1)*x(i) + x(i)*x(i)*y(i);
48 dydt(i) = (y(i-1)-2*y(i)+y(i+1))*cy + B*x(i) - x(i)*x(i)*y(i);
49 end
50
51 dfdt = [dxdt; dydt];
52
53 % --------------------------------------------------------------------------
54
55 function [tspan,y0,options] = init(N)
56 tspan = [0; 10];
57 A = 2;
118
58 B = 4.6;
59
60 y0 = zeros(2*N,1);
61
62 for i=1:N
63 y0(i) = A + 2*sin(pi*i/(N+1));
64 y0(N+i) = B/A - 0.5*sin(pi*i/(N+1));
65 end
66 handles = feval(@bruss);
67 options = odeset(’Vectorized’,’on’, ’Jacobian’, handles(3), ’JacobianP’, handles(4));
68
69 % --------------------------------------------------------------------------
70
71 function dfdxy = jacobian(t,y,N,L)
72 x = y(1:N);
73 y = y(N+1:2*N);
74 A = 2;
75 B = 4.6;
76 Dx = 0.0016;
77 Dy = 0.008;
78 x0 = A; x1 = A;
79 y0 = B/A; y1 = B/A;
80 L2 = L^2;
81 h = 1/(N+1);
82 cx = (Dx/L2)/(h*h);
83 cy = (Dy/L2)/(h*h);
84
85
86 %
87 % Sparse jacobian
88 %
89 A=zeros(2*N,3);
90 A(1:N-1,2)=cx;
91 A(1:N,3)=-2*cx -(B+1) + 2*x(1:N).*y(1:N);
92 A(1:N,4)=cx;
93
94 A(N+1:2*N,2) = cy;
95 A(N+1:2*N,3) = -2*cy -x(:).*x(:);
96 A(N+2:2*N,4) = cy;
97
98
99 A(1:N,1) = B - 2*x(:).*y(:);
100 A(N+1:2*N,5) = x(:).*x(:);
101
102 dfdxy = spdiags(A, [-N,-1:1,N] , 2*N, 2*N);
103 return
104
105 %
106 % Full matrix
107 %
108 dfxdx = zeros(N,N);
109 dfydy = zeros(N,N);
110 dfxdy = zeros(N,N);
111 dfydx = zeros(N,N);
112
113 for i=1:N
114 if i>1, dfxdx(i,i-1) = cx; end;
119
115 dfxdx(i,i) = -2*cx -(B+1) + 2*x(i)*y(i);
116 if i<N, dfxdx(i,i+1) = cx; end;
117 dfxdy(i,i) = x(i)*x(i);
118
119 if i>1, dfydy(i,i-1) = cy; end;
120 dfydy(i,i) = -2*cy -x(i)*x(i);
121 if i<N, dfydy(i,i+1) = cy; end;
122 dfydx(i,i) = B - 2*x(i)*y(i);
123 end
124
125 dfdxy = [ dfxdx, dfxdy; dfydx, dfydy ];
126
127 % --------------------------------------------------------------------------
128
129 function dfdp = jacobianp(t,y,N,L)
130 x = y(1:N);
131 y = y(N+1:2*N);
132 A = 2;
133 B = 4.6;
134 Dx = 0.0016;
135 Dy = 0.008;
136 x0 = A; x1 = A;
137 y0 = B/A; y1 = B/A;
138 L2 = L^2;
139 h = 1/(N+1);
140 cx = (Dx/L2)/(h*h);
141 cy = (Dy/L2)/(h*h);
142 kx = (-2/L)*cx;
143 ky = (-2/L)*cy;
144
145 Sx = zeros(N,1);
146 Sy = zeros(N,1);
147
148 Sx(1) = kx*(x0-2*x(1)+x(2));
149 Sy(1) = ky*(y0-2*y(1)+y(2));
150
151 Sx(N) = kx*(x(N-1)-2*x(N)+x1);
152 Sy(N) = ky*(y(N-1)-2*y(N)+y1);
153
154 i=2:N-1;
155 Sx(i) = kx*(x(i-1)-2*x(i)+x(i+1));
156 Sy(i) = ky*(y(i-1)-2*y(i)+y(i+1));
157
158 dfdp = [ zeros(2*N,1) [Sx;Sy] ];
159
160 % --------------------------------------------------------------------------
bruss.m
The model is implemented with 2 parameters: N and L; the values of A, B, Dx , Dy are hard-
coded. Note that N is a parameter that cannot vary during the continuation. Therefore it
does not have entries in Jacobianp. We should let the pde 1 curve know that bruss.m is the
active file, the initial values of the parameters N and L are respectively 20 and 0.06 and the
active parameter is L, i.e. the second parameter of bruss.m. So, first of all we have to get
the approximate equilibrium solution which is provided in a subfunction ’init’ of ’bruss.m’.
120
2.8
2.6
LP
2.4
x(20)
2.2
LP
1.8
1.6
0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26
L
Figure 33: An quilibrium curve of bruss.m
We first locate the handle to the subfunction init and call it.
>N = 20; L = 0.06;

>handles = feval(@bruss);
>[t,x0,options] = feval(handles{1},N);
It sets the number of state variables to 2N and makes an initial vector x0 of length 2N
containing the values of the approximate equilibrium solution. Now we inform the pde 1
curve that the second parameter of bruss.m is the active parameter and what the default
values of the other parameters are. We also set some options.
>[x1,v1] = init_EP_EP(@bruss,x0,[N;L], [2]);

>opt = contset;opt=contset(opt,’MinStepsize’, 1e-5);
>opt=contset(opt,’MaxCorrIters’, 10);
>opt=contset(opt,’MaxNewtonIters’, 20);
>opt=contset(opt,’FunTolerance’, 1e-3);
>opt=contset(opt,’Singularities’,1);
>opt=contset(opt,’MaxNumPoints’,500);
>opt=contset(opt,’Locators’,[]);
We start the continuation process by the statement

[x,v,s,h] = cont(@pde 1,x1,v1,opt).
In this case the number of state variables can be a parameter and the Jacobian can be sparse.
Using the command
cpl(x,v,s,[41;20]);
we can plot the 20th component of x as a function of L (which is the 41th continuation
variable), see Figure 33.
The above test can be executed by running testbrusselator; this file is in the directory
Testruns.
121
References
[1] Bashir Al-Hdaibat, Willy Govaerts, Yuri A. Kuznetsov and Hil G.E.
Meijer, Initialization of homoclinic solutions near Bogdanov-Takens points:
Lindstedt-Poincaré compared with regular perturbation method, SIAM J. Appl.
Dyn. Syst. 15(2). p.952-980 (2016).
[2] E.L. Allgower and K. Georg, Numerical Continuation Methods: An introduc-
tion, Springer-Verlag, 1990 .
[3] W.J. Beyn, A. Champneys, E. Doedel, W. Govaerts, Yu.A. Kuznetsov,
and B. Sandstede, Numerical continuation and computation of normal forms. In:
B. Fiedler, G. Iooss, and N. Kopell (eds.) “Handbook of Dynamical Systems : Vol
2”, Elsevier 2002, pp 149 - 219.
[4] Champneys, A.R. and Kuznetsov Yu.A. 1994. Numerical detection and con-
tinuation of codimension-two homoclinic orbits. Int. J. Bifurcation Chaos, 4(4),
785-822.
[5] Champneys, A.R., Kuznetsov Yu.A. and Sandstede B. 1996. A numerical
toolbox for homoclinic bifurcation analysis. Int. J. Bifurcation Chaos, 6(5), 867-887.
[6] C. De Boor and B. Swartz, Collocation at Gaussian points, SIAM Journal on
Numerical Analysis 10 (1973), pp. 582-606.
[7] Demmel, J.W., Dieci, L. and Friedman, M.J. 2001. Computing connecting
orbits via an improved algorithm for continuing invariant subspaces. SIAM J. Sci.
Comput., 22(1), 81-94.
[8] V. De Witte, W. Govaerts, Yu. A. Kuznetsov and M. Friedman, In-
teractive Initialization and Continuation of Homoclinic and Heteroclinic Orbits in
MATLAB, ACM Transactions on Mathematical Software. Volume 38, Issue 3, Ar-
ticle Number: 18, DOI: 10.1145/2168773.2168776 Published: APR 2012
[9] V. De Witte, F. Della Rossa, W.Govaerts and Yu.A. Kuznetsov, Nu-
merical Periodic Normalization for Codim2 Bifurcations of Limit Cycles: Compu-
tational Formulas, Numerical Implementation, and Examples, SIAM J. Applied
Dynamical Systems 12,2 (2013) 722-788. DOI: 10.1137/120874904
[10] A. Dhooge, W. Govaerts and Yu. A. Kuznetsov, MATCONT : A matlab
package for numerical bifurcation analysis of ODEs, ACM Transactions on Mathe-
matical Software 29(2) (2003), pp. 141-164.
[11] A. Dhooge, W.Govaerts, Yu. A. Kuznetsov, H. G. E. Meijer and B.
Sautois: New features of the software MatCont for bifurcation analysis of dynam-
ical systems, Mathematical and Computer Modelling of Dynamical Systems, Vol.
14(2), pp. 147-175, Published: 2008.
https://doi.org/10.1080/13873950701742754.
[12] E. Doedel and J Kernévez, AUTO: Software for continuation problems in or-
dinary differential equations with applications, California Institute of Technology,
Applied Mathematics, 1986.
122
[13] Doedel, E.J. and Friedman, M.J.: Numerical computation of heteroclinic or-
bits, J. Comp. Appl. Math. 26 (1989) 155-170.
[14] E.J. Doedel, A.R. Champneys, T.F. Fairgrieve, Yu.A. Kuznetsov, B.

Sandstede and X.J. Wang, auto97-00 : Continuation and Bifurcation Soft-
ware for Ordinary Differential Equations (with HomCont), User’s Guide, Concordia
University, Montreal, Canada (1997-2000). (http://indy.cs.concordia.ca).
[15] Doedel, E.J., Govaerts W., Kuznetsov, Yu.A.: Computation of Periodic So-
lution Bifurcations in ODEs using Bordered Systems, SIAM Journal on Numerical
Analysis 41,2(2003) 401-435.
[16] Doedel, E.J., Govaerts, W., Kuznetsov, Yu.A., Dhooge, A.: Numerical
continuation of branch points of equilibria and periodic orbits, Int. J. Bifurcation
and Chaos, 15(3) (2005), 841-860.
[17] Ermentrout, B.: Simulating, Analyzing, and Animating Dynamical Systems.

Siam Publications, Philadelphia, 2002.
[18] Freire, E., Rodriguez-Luis, A., Gamero E. and Ponce, E., A case study for
homoclinic chaos in an autonomous electronic circuit: A trip form Takens-Bogdanov
to Hopf- Shilnikov, Physica D 62 (1993) 230–253.
[19] Friedman, M., Govaerts, W., Kuznetsov, Yu.A. and Sautois, B. 2005.
Continuation of homoclinic orbits in matlab. Lecture Notes in Computer Science,
3514, 263-270.
[20] Genesio, R. and Tesi, A. Harmonic balance methods for the analysis of chaotic
dynamics in nonlinear systems. Automatica 28 (1992), 531-548.
[21] Genesio, R., Tesi, A., and Villoresi, F., Models of complex dynamics in
nonlinear systems. Systems Control Lett. 25 (1995), 185-192.
[22] W.J.F. Govaerts, Numerical Methods for Bifurcations of Dynamical Equilibria,

SIAM, 2000.
[23] Govaerts, W. and Sautois, B.: Phase response curves, delays and synchroniza-
tion in matlab. Lecture Notes in Computer Science, 3992 (2006), 391-398.
[24] Govaerts, W. and Sautois, B.: Computation of the phase response curve: a
direct numerical approach. Neural Comput. 18(4) (2006), 817-847.
[25] Yu. A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag,

1998. (third edition 2004).
[26] Yu. A. Kuznetsov and V.V. Levitin, CONTENT: Integrated En-

vironment for analysis of dynamical systems. CWI, Amsterdam 1997:
ftp://ftp.cwi.nl/pub/CONTENT
[27] matlab, The Mathworks Inc., http://www.mathworks.com.
123
[28] Yu. A. Kuznetsov, W. Govaerts, E.J. Doedel and A. Dhooge, Numerical
periodic normalization for codim 1 bifurcations of limit cycles, SIAM J. Numer.
Anal. 43 (2005) 1407-1435.
[29] Yu.A. Kuznetsov, H.G.E. Meijer, W. Govaerts and B. Sautois, Switching

to nonhyperbolic cycles from codim 2 bifurcations of equilibria in ODEs, Physica
D 237 No. 23 (2008) 3061-3068 (ISSN 0167-2789).
[30] Kuznetsov, Yu.A., Meijer H.G.E., Al Hdaibat, B. and Govaerts, W., Im-
proved homoclinic predictor for Bogdanov-Takens Bifurcation, International Jour-
nal of Bifurcations and Chaos, 24(4) (2014). Article Number: 1450057. DOI:
10.1142/S0218127414500576
[31] Yu.A. Kuznetsov, H.G.E. Meijer, B. Al-Hdaibat and W. Govaerts, Ac-

curate Approximation of Homoclinic Solutions in Gray-Scott Kinetic Model. In-
ternational Journal of Bifurcation and Chaos, Volume: 25(9) August 2015. Article
Number: 1550125
DOI: 10.1142/S0218127415501254
[32] W. Mestrom, Continuation of limit cycles in matlab, Master Thesis, Mathemat-

ical Institute, Utrecht University, The Netherlands, 2002.
[33] Morris, C., Lecar,H., Voltage oscillations in the barnacle giant muscle
fiber,Biophys J. 35 (1981) 193–213.
[34] N. Neirynck, Advances in numerical bifurcation software: MatCont. PhD thesis,

Ghent University, Belgium 2019. https://biblio.ugent.be/publication/8615817.
[35] A. Riet, A Continuation Toolbox in matlab, Master Thesis, Mathematical Insti-

tute, Utrecht University, The Netherlands, 2000.
[36] C. Stéphanos, Sur une extension du calcul des substitutions linéaires, J. Math.
Pures Appl. 6 (1900) 73-128.
[37] Terman, D., Chaotic spikes arising from a model of bursting in excitable mem-
branes, Siam J. Appl. Math. 51 (1991) 1418–1450.
[38] Terman, D., The transition from bursting to continuous spiking in excitable mem-
brane models, J. Nonlinear Sci. 2, (1992) 135–182.
124

Manual Mat Con Taug 2019

Uploaded by

Copyright:

Available Formats

Manual Mat Con Taug 2019

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Manual Mat Con Taug 2019

Uploaded by

Copyright:

Available Formats

MATCONT:

Continuation toolbox for ODEs in

August 2019, adapted for version MatCont 7.1.

2 Mathematical aspects of numerical continuation and handling of singulari-

3 General software aspects of MatCont 18

5 Time integration and Poincaré maps 41

7 Continuation of limit cycles 53

8 Continuation of codim 1 bifurcations 68

9 Continuation of codim 2 bifurcations 98

10 Continuation of homoclinic and heteroclinic orbits 107

A Continuer example: a curve object 115

B The Brusselator example: Continuation of a solution to a boundary value

• Allow for easy extendibility.

2 BP CP BT ZH HH GH CPC BPC R1 R3 R4 CH LPNS PDNS R2 NSNS LPPD GPD

An arrow in Figure 1 from O to EP or LC means that by starting time integration from

2 NSS NSF NFF DR* ND* TL* SH OF* IF* NCH

1.2 Testruns and how to start with Cl MatCont

• testodefile.m, discussed in §4.3 (basic use of the system definition file)

• testtimeinteg.m and testtimeintegJacobian, discussed in §5.1.3 (time integration)

• testPoincare.m, discussed in §5.2.1 (Poincaré map)

• testadaptPRC.m, discussed in §7.8 (continuation of equilibria, detection of a Hopf point,

• testselectcycle.m, discussed in §7.4 (starting the continuation of limit cycles from

• testbratu.m and testbratu2.m, discussed in §6.5 (user function, equilibrium contin-

• testequilcataloscill.m and testLPcataloscill.m, discussed in §8.1.5 (equilibrium

• testLPHopfcataloscill.m, discussed in §8.2.5 (equilibrium continuation, detection of

neutral equilibria, detection of GH points)

• testadapt.m, discussed in §7.7 (continuation of equilibria, detection of Hopf points)

• testadapt2.m, discussed in §7.7 (starting a period doubled orbit from a PD point;

• testadapt3.m, discussed in §8.3.5 (starting a curve of period doubling bifurcations from

• testEquilMLfast.m, discussed in §8.4.4 (continuation of equilibria, detection of limit

• testLCMLfast.m, discussed in §8.4.4 (continuation of limit cycles and detection of an

• testLPCMLfast.m, discussed in §8.4.4 (continuation of an LPC curve with two free

• testtorBPC1.m and testtorBPC2.m, discussed in §8.5.4 (detection of a Neimark-Sacker

• testtorBPC3.m, discussed in §8.5.4 (continuation of Neimark-Sacker bifurcations in two

• testtorBPC4.m, discussed in §9.2.4 (continuation of equilibria, detection of a Hopf point

• testtorBPC5.m, discussed in §9.2.4 (continuation of limit cycles, detection of one LPC

• testtorBPC6.m, discussed in §9.2.4 (branching off in a BPC, detection of a NS and a

• testtorBPC7.m, discussed in §9.2.4 (continuation of a BPC with three free parameters)

• testmyml.m and homoc1.m, discussed in §10.5 (computation of a curve of equilibria

• testdrawcurve.m, discussed in Appendix A (drawing a curve by using a continuation

• testbrusselator.m, discussed in Appendix B (continuation of an equilibrium solution

1.3 Computational routines in MatCont

Computation of the normal form coefficients of bifurcations is another important task.

1.5 Software requirements

• correction of the predicted point

The question is how to choose the function g(x).

2.2.1 Pseudo-arclength continuation

X k+1 = X k − Hx−1 (X k )H(X k ) (4)

2.2.2 Moore-Penrose continuation

Definition 1 Let A be an N × (N + 1) matrix with maximal rank. Then the Moore-Penrose

min{||x − X 0 || | F (x) = 0} (11)

2.3 Stepsize control

2.4.1 Test functions

F (x) = 0, F : Rn+1 → Rn . (23)

The singularity S will then be detected if

φ(xi )φ(xi+1 ) < 0 . (24)

2.4.2 Multiple test functions

∀j ∈ [1, nt ] : φj (xi )φj (xi+1 ) < 0 (27)

∀j ∈ [1, nt ] : x∗ = x∗j and φj (x∗j ) = 0 (28)

2.4.4 User location