Continuation and Bifurcation in Nonlinear PDEs – Algorithms, Applications, and Experiments

Numerical continuation and bifurcation methods can be used to explore the set of steady and time–periodic solutions of parameter dependent nonlinear ODEs or PDEs. For PDEs, a basic idea is to first convert the PDE into a system of algebraic equations or ODEs via a spatial discretization. However, the large class of possible PDE bifurcation problems makes developing a general and user–friendly software a challenge, and the often needed large number of degrees of freedom, and the typically large set of solutions, often require adapted methods. Here we review some of these methods, and illustrate the approach by application of the package pde2path to some advanced pattern formation problems, including the interaction of Hopf and Turing modes, patterns on disks, and an experimental setting of dead core pattern formation.


Introduction
Ordinary differential equation (ODE) models and partial differential equation (PDE) models usually come with a number of parameters λ (growth factors, diffusion constants . . . ), and an important task is to characterize the dependence of solutions on parameters. Numerical continuation aims to trace branches of solutions u of nonlinear equations G(u, λ) = 0 ∈ Y through parameter space, where X, Y are Banach spaces and G : X × R → Y has some smoothness detailed below, and where s is an arclength along the solution branch, which may fold back and forth in λ. If u is a steady solution (fixed point) of an ODE d dt u = −G(u, λ), then X = Y = R n and G(u, λ) = 0 is an algebraic system, and if u : → R N is a steady solution of a PDE ∂ t u = −G(u, λ), ⊂ R d a bounded domain, then X and Y are function spaces and G(u, λ) = 0 is a boundary value problem. Additionally, often H. Uecker hannes.uecker@uol.de 1 Institut für Mathematik, Carl von Ossietzky Universität Oldenburg, Oldenburg, Germany time periodic orbits u = u(t) with u(t + T ) = u(t) of ODEs and PDEs are of special interest.
If ∂ u G(u(s 0 ), λ(s 0 )) ∈ L(X, Y ) is invertible, then the implicit function theorem yields that the branch (s 0 − δ, s 0 + δ) s → (u(s), λ(s)) is locally unique. Of special interest are singular points where ∂ u G(u(s 0 ), λ(s 0 )) is not invertible and hence other branches may bifurcate. At such points we may aim to switch to bifurcating branches, obtaining so-called bifurcation diagrams, altogether aiming at an as complete as possible description of the set of solutions of a given problem.
Numerical continuation and bifurcation methods have been an active field over the last 50 years, see, e.g., [32,37,45,87] for textbooks and research monographs, and, e.g., [23,24,52,60,65] for reviews. One main starting point is the (pseudo-)arclength parametrization of solution branches [36], leading to algorithms implemented in many packages such as AUTO [25], matcont [19], xppaut [27], aimed (primarily) at the numerical bifurcation analysis of ODEs, and coco [14], which is a rather general toolbox. See also [79] for a recent Julia package. For PDEs, the first step is often a spatial discretization, which turns the PDE into a large ODE system, such that in principle the above software can be applied. However, the possibly large number of degrees of freedom (DoF), in particular for spatially two-or three-dimensional problems, often requires adapted algorithms. Software aimed at PDEs includes LOCA [59] and oomph [33].
Here we give an overview of the numerical computation of bifurcation diagrams for nonlinear PDEs, first reviewing some basics, and then giving examples using the MATLAB package pde2path [68,73]. This is aimed at PDEs of the form where M d ∈ R N ×N is the (dynamical) mass matrix, M d = Id in many cases, and with u=u(x, t)∈R N (N components), x∈ ⊂R d some bounded domain, d∈{1, 2, 3} (the 1D, 2D and 3D case, respectively), and time t ≥ 0, and where (2) can be complemented with various boundary conditions (BCs). In (3), c is a diffusion tensor, b describes advection, and a ∈ R N ×N and f ∈ R N should be thought as describing linear and nonlinear terms without spatial derivatives. The mass matrix M d ∈ R N ×N in (2) may be singular, and this gives much flexibility to (2), for instance to also implement parabolic-elliptic coupled problems. All the tensors/vectors in (3) can depend on x and u, although mostly we restrict to semilinear problems, where c does not depend on u. A focus of pde2path is on bifurcations of solution branches of the steady state problem but we also consider Hopf bifurcations and time-periodic orbits and their bifurcations in (2). The default setting of pde2path uses the finite element method (FEM), essentially as provided by the package OOPDE [58], to first discretize u in space, in 1D, 2D, or 3D, using Lagrangian P 1 elements. However, some higher order FEM is also provided, and for instance [71] contains examples of FEM-free implementation of various right hand sides G, mostly based on Chebychev and Fourier methods. The three main design goals of pde2path are versatility, easy usage, and easy customization. Therefore, in [67,76] and [68] we explain the usage of pde2path via a rather large number of working examples included in the software download at [73], where we also provide tutorials which include further comments on implementations.
Here we proceed similarly, in a condensed form: In §2 we explain some basic notions of bifurcations, and a few algorithms of numerical continuation and bifurcation. In §3 we then give some examples of applications of pde2path, first to rather standard PDE bifurcation problems, which however do not seem to have been treated in this way before, namely: 1. In §3.1 we consider the Reaction-Diffusion type system over the interval = (−l, l) with Neumann BCs for u 1 and u 2 , taken from [53].
The system describes charge transport in a layered semiconductor; u 1 is an interface charge in the device, u 2 a normalized voltage across it, and (τ, j 0 , α, d) ∈ R 4 is a parameter vector. There is the trivial (spatially homogeneous) branch , u * 2 = j 0 + u * 1 , which may undergo Turing instabilities (bifurcation to finite wave number spatially periodic steady states) and Hopf instabilities (bifurcation to spatially homogeneous time-periodic orbits). In particular, there are codimension-two points, near which both instabilities occur, and their interaction may lead to a variety of states, including Turing-Hopf mixed modes, i.e. localized Hopf modes embedded in Turing structures, or vice versa. Similar results have been obtained in [15] for the so-called Brusselator system instead of (5). Some of these patterns can be studied via amplitude equations, but for larger amplitudes [53] and [15] use direct numerical simulation (DNS), aka time integration. Here we shall study these patterns via numerical continuation and bifurcation in a complementing way. 2. Following [81], we consider the cubic-quintic Swift-Hohenberg (SH) equation on a disk domain = {(x 1 , x 2 ) : x = x 2 1 + x 2 2 < R}, where = ∂ 2 x 1 + ∂ 2 x 2 is the Laplacian, with Neumann BCs ∂ n u = ∂ n u = 0, where n is the outer normal at x = R. Equations of type (6) (with various nonlinearities, for instance f (u) = νu 2 − u 3 instead of f (u) = νu 3 − u 5 in (6)) are prototypical examples for finite wave number (Turing) pattern formation, and are thus also studied as model problems in [68], over 1D, 2D and 3D boxes (intervals, rectangles and cuboids). The new feature in [81] is the disk domain. The problem has the trivial branch u ≡ 0, and for small ε > 0 we find various subcritical bifurcating patterns, and also study some of their secondary (and tertiary) bifurcations.
with f (u) ∼ u γ + at u = 0, u + = max(u, 0), where 0 < γ < 1 and hence f (u) is not differentiable (and not even Lipschitz) at u = 0. For large λ, such equations can have dead core solutions which feature a subdomain 0 ⊂ where u ≡ 0, which is not possible for Lipschitz nonlinearities. The non-differentiability of f means that standard bifurcation analysis does not apply to (7). However, using minor modifications of the standard pde2path setup we can still treat (7) and generalizations of (7) numerically, and find various branches with dead core solutions.
Our main aim with these examples is to present interesting bifurcation problems and how they can be studied by numerical continuation and bifurcation. Other recent applications of pde2path are given in, e.g., [6,9,26,41,44,55,62,69,75,78,82,83,86]. In Table 1 we collect some acronyms frequently used in PDEs and bifurcation theory, and some specific for this paper. Initially, pde2path was planned as an in-house tool to quickly enable students to study steady state bifurcations for 2D reaction-diffusion systems, based on the MAT-LAB pdetoolbox (see also [73] for octave compatibility). As it was well received by students and colleagues, in 2013 we decided to go public with a basic version, after just a few months of development, and only treating 2D steady state problems (and simple bifurcation points). At the time I was not quite aware of what we were getting into. Ever since, I have learnt a lot about bifurcations, the FEM, and software organization, and one experience I'd like to share is: even if you just have a small starting version of a software that might be useful, try to make it simple for others, provide working examples, and go public. If the software turns out to be useful, the feedback and requests by users are a tremendous help and motivation. Therefore, also thanks to all the users, and please keep sending bug reports, feedback, and requests.

Continuation and Bifurcation
In [68, Part I], I give a brief discussion of the theoretical background underlying continuation and bifurcation. For convenience, here I summarize the main ideas; although none of this is original, I give only a few references, and instead refer to [68] for further background and references. Readers familiar with (numerical) bifurcation theory can skip directly to §3 and come back to §2 when needed.
The fundamental result underlying all bifurcation theory is the Implicit Function Theorem (IFT), which can be proved via the contraction mapping theorem. See, e.g., [85, §4.7] for precise statements and a thorough discussion, while here we just recall the main ideas. As an example, consider For each λ ≥ −1 we have the two solutions u = u ± (λ) = 1 ± √ 1 + λ, and together we have the solution branch showing that here it is useful to see λ as a function of u. However, all we want to use now is that for λ 0 = 0 we can "guess" the solution u 0 = 0. We treat λ as a parameter and seek zeros of G(u, λ) near (u, λ) = (u 0 , λ 0 ) = (0, 0) in form of a (smooth) resolution u = u(λ). Implicit differentiation yields If ∂ u G(u 0 , λ 0 ) = 0, then we can solve for and altogether we can formally obtain the Taylor expansion of u(λ) around λ 0 . In the example we have u (0) = −1/2, u (0) = 1/4, and hence Two natural question are: a) Does the Taylor expansion (12) converge (near λ = 0)? Are there other solutions of G(u, λ) = 0 near (u, λ) = (0, 0), not given by u(λ)? Of course, for our example we know all the answers because we explicitly know the solution: Clearly, (12) can only converge locally, as near (u, λ) = (1, −1) there does not exists a resolution u = u(λ). Instead, a so-called saddle-node bifurcation occurs, which is also called a fold or a turning point. The theory behind (12) (and much more) is summarized in the IFT: Let X, Y and be Banach spaces, and let B X ε (u 0 ) denote the open ε-ball around u 0 in X, and L(X, Y ) the Banach space of continuous linear operators A : X → Y .
. Then the following holds:

Remark 2.2 a)
The results in ii) and iii) justify Taylor expansions, and the proof of Theorem 2.1 is constructive in the sense that it yields formulas for the Taylor expansion of H , like in (12). b) Concerning example (8), i) gives the local uniqueness, and ii) gives the local convergence. From the proof of the IFT (or directly from the calculus leading to (12)) we can estimate the convergence radius, namely |λ| < 1, but in practice we are satisfied with the local results. In particular, the branch (u(λ), λ) through (0, 0) can be continued as long as ∂ u G(u, λ) = 0. As already indicated in (9), it may be useful to switch the roles of u and λ. In (8) we have ∂ λ G(u, λ)= − 1 =0 (independent of (u, λ)), and hence we can always find a (local, and in fact global) resolution λ=λ(u) (namely (9)).

Standard Bifurcation Examples for ODEs
To explain the notion of bifurcation and the meaning of bifurcation diagrams, we consider the three so-called normal forms of elementary steady bifurcations, and a simple Hopf bifurcation example.

Example 2.3 a)
Saddle-node (or fold) bifurcation, G(u, λ) = λ − u 2 . Here the two solution branches u = u ± (λ) = ± √ λ exist for λ ≥ 0, see Fig. 1(a1). Of course, equivalently λ = u 2 , but for now we stick to the u = u(λ) point of view. For any λ > 0 we have ∂ u G(u ± (λ), λ) = ±2 √ λ = 0, such that the solutions are locally unique by the IFT. Plots of solutions (or of some functionals of solutions, e.g., norms) as functions of the parameter as in Fig. 1(a1) are called bifurcation diagrams (BDs). In (a2) we show the same branches and additionally sketch the flow of the associated ODE d dt u = G(u, λ). Moreover, as custom in bifurcation diagrams, stable branches (see below) are plotted in full (or thick) lines, while unstable branches as dashed (or thin) lines.
The point (u, λ) * = (0, 0) in all three cases is called a bifurcation point. In b) and c) the point (u, λ) * is also called a branch point, as two branches intersect there. At the saddle-node bifurcation there is no "branching off", and we rather prefer to call (u, λ) * a fold point, or simply a fold. Importantly, the transcritical bifurcation and the pitchfork bifurcation are non-generic. This means that they require some special structure (or symmetries) of a nonlinear problem G(u, λ) = 0 to occur at all, but such symmetries are often enforced by physics. See [68, §1.2] (and the references therein) for more on genericity (or structural stability), co-dimensions, and imperfections.
By the IFT, a necessary condition for all three bifurcations is that ∂ u G(u * , λ * ) is not invertible, i.e., has a zero eigenvalue. Additionally, the eigenvalues of the Jacobian A = ∂ u G(u * , λ * ) at a steady state u * often determine the stability of u * . The fixed point u * is called stable for d dt u = −G(u, λ), if solutions to nearby initial conditions u 0 stay close to u * , i.e., if for all ε > 0 there exists a δ > 0 such that u 0 − u * < δ implies u(t) − u * < ε. If additionally u(t) − u * → 0, then u * is called asymptotically stable. If u * is not stable, then it is called unstable. Bifurcating branches are often classified as super-or subcritical. Supercritical means that the bifurcating branch exists in the λ-range where the original (trivial) branch has lost stability, while subcritical means that it exists where the trivial branch is stable.
If u * is a fixed point of d dt u = −G(u, λ), then we may consider the linearization (13) at u * . This can be solved explicitly by an exponential ansatz v(t) = ce −μt φ, where μ ∈ C is an eigenvalue of A and φ the associated eigenvector (for semisimple eigenvalues, i.e., same algebraic and geometric multiplicity, to be generalized in case that A has Jordan blocks), and from this it follows that u * is (i) asymptotically stable, if all eigenvalues μ of A fulfill Reμ > 0; (ii) unstable, if at least one eigenvalue μ has Reμ < 0, or if there are Jordan blocks to a zero eigenvalue.
The remaining cases, that there is a semisimple eigenvalue μ = 0, or that there are purely imaginary eigenvalues μ = ±iω, are associated to possible bifurcations, and the dynamics close to u * must be studied in detail, for instance via center-manifold reduction. A similar principle of linearized stability also holds for many evolutionary PDEs, see, e.g., [61,Theorem 5.2.23].
For λ > 0 we obtain the (stable) fixed point r = √ λ of the d dt r equation, and hence a family of (stable) periodic orbits (POs) bifurcates from the trivial branch u=0 at λ=0, see Fig. 2. This is called a Hopf bifurcation. For fixed λ > 0 the family attracts every solution with an exponential rate O(exp(−2λt)).

Remark 2.5
The stability of POs is somewhat less straightforward to define (via Poincaré sections and maps) than the stability of steady states, and more difficult to analyze (via Floquet multipliers). Similarly, bifurcations from POs (nontrivial Floquet multipliers μ with |μ| = 1) are more difficult, both analytically and numerically, than steady bifurcations, or (Hopf-) bifurcations of POs from steady states, and we refer to [68, Chap. 1] for details and references.  (14). In the BD we now need to plot some norm of u (Color figure online)

The Crandall-Rabinowitz Theorem, and Remarks
From the examples in §2.1 we may expect that the crossing of a simple eigenvalue through 0 generically leads to steady bifurcations, while a pair of nonzero eigenvalues crossing the imaginary axis leads to Hopf bifurcations. This is true, under some technical assumptions, and can be proved via Liapunov-Schmidt reduction, where again we refer to [68, §2.2], and here only state a prototype bifurcation theorem for steady bifurcation from simple eigenvalues [11]. 1 We consider (4), i.e., G(u, λ) = 0, and for simplicity assume that this has the trivial branch u = 0, λ ∈ R, and without loss of generality consider a bifurcation at (u, λ) Theorem 2.6 (Crandall-Rabinowitz.) Consider (4) and assume that additional to (H1),(H2) we have the transversality condition where (4), and the nontrivial branch bifurcating at (0, 0) is given locally by a C 1 curve where u(s) = sφ + O(s 2 ). All solutions of G(u, λ) = 0 in a neighborhood of (0, 0) are either on γ or on the trivial branch.
Like the proof of the IFT, the proof of Theorem 2.6 is constructive in the sense that it yields bifurcation formulas, i.e., formulas for λ (0), λ (0), . . . (if desired). Moreover, there are various extensions of Theorem 2.6, concerning, e.g., the exchange of stability at a BP, and possible global behaviors of bifurcating branches (Rabinowitz alternative).
The assumption that 0 is a simple eigenvalue is crucial; in case of higher multiplicity, different things can happen: no bifurcation, or bifurcation of higher dimensional 'leaves', or bifurcation of several branches (often more than the multiplicity of the zero eigenvalue). In principle, one must then derive and solve the so-called (leading order) algebraic bifurcation equations. In this, often symmetry (which usually creates the higher multiplicity of the eigenvalue) can be heavily exploited to simplify the problem. This comes under the name of equivariant branching, and is often crucial for bifurcations in PDEs, but here we refer to [31,34] and the references therein for details, and [68, §2.5] for an overlook.

Basic Numerical Algorithms
The general bifurcation theory, of which Theorem 2.6 is one example, works in Banach spaces, while the examples in §2.1 were for simple scalar (resp. 2-component for the Hopf case) ODEs. For the numerical continuation and bifurcation in pde2path we first discretize PDEs (of the form (2)) in space, yielding (high dimensional) algebraic problems (with a slight abuse of notation) or ODEs (or, if the dynamical mass matrix M d ∈ R n×n is singular, differential algebraic systems) For the computation of branches (and bifurcations) for (17) and (18), some special parametrizations and methods have been established. Here we recall the basic ideas of arclength continuation, focusing on methods which are implemented in pde2path.
Besides the algorithms for (one parameter) branch continuation and steady bifurcations, we shall give some remarks on multiple-parameter continuation of, e.g., fold points (FPs), branch points (BPs), and Hopf bifurcation points (HPs). To some extent, these methods can be formulated in general Banach spaces X, but for convenience we restrict to X = R n . In particular, ∂ u G(u, λ) ∈ R n×n is always Fredholm of index 0. We essentially present the algorithms in an abstract mathematical way, but where useful add comments pertaining to the pde2path implementation.
Remark 2.7 a) We emphasize that the notation G(u, λ) = 0 ∈ R n is motivated by u ∈ R n corresponding to the values of the field u : → R N at the n p mesh points, hence n = Nn p , while λ ∈ R is one "active" scalar parameter. Strictly speaking, for the computation of solution branches s → (u(s), λ(s)) ∈ R n+1 there is no reason to distinguish (in notation or in substance) between the field values u = (u 1 , . . . , u n ) ∈ R n and the parameter λ ∈ R. The only thing that matters is that we have n equations in the n + 1 unknowns (u, λ). Thus, we could as well rename λ = u n+1 , and essentially this is what is done in pde2path, i.e., all parameters are appended to the vector u = (u, pars), and the active parameter "λ" (or several active ones) is (are) marked by a pointer. b) Often we meet situations where the system of PDEs must be extended by additional equations, for instance so-called phase conditions (PCs). Each additional such equation Q i = 0, i = 1, . . . , q, naturally requires to free ("activate") one more parameter, such that we obtain one primary active parameter λ, and a number of secondary active parameters w ∈ R q . Similarly, the continuation of POs typically requires to free at least one additional active parameter, usually the period T , and to set a PC as at least one additional equation. In principle, we could put the additional equations into the present framework by extending u to (u, λ, w) ∈ R n+1+q and G to (G, Q 1 , . . . , Q q ) : R n+1+q × R → R n+q , but in practice we find it more transparent to keep the distinction between the field u and the parameter(s) (λ, w), and between the "PDE" G(u, λ) = 0 or ∂ t u = −G(u, λ) and the auxiliary equations Q 1 = 0, . . . , Q n q = 0. Thus, we essentially first focus on the 1-parameter case (q = 0), and come back to fold-point-, branch-point-, and Hopf-point-continuation in §2.4.

Arclength Continuation A standard method for numerical continuation of branches
parameterized by s ∈ R and the extended system where p is used to make s an approximation of arclength on the solution branch. Given s 0 and a point (u 0 , λ 0 ) := (u(s 0 ), λ(s 0 )), and additionally a tangent vector Here 0 < ξ < 1 is a weight, typically chosen as ξ = 1/n u , and τ 0 is assumed to be normalized in the weighted norm For fixed s and τ 0 ξ = 1, p(u, λ, s) = 0 thus defines a hyperplane perpendicular (in the inner product ·, · ξ ) to τ 0 at distance ds := s − s 0 from (u 0 , λ 0 ). We may then use a predictor The idea of arclength continuation for a solution of (19) on that hyperplane, followed by a corrector using Newton's method in the form where z = A −1 b stands for the solution of the linear system Az = b, see Fig. 3. 2 Since ∂ s p = −1, on a smooth solution arc we have Thus, after convergence of (22b) yields a new point (u 1 , λ 1 ) with Jacobian A, the tangent direction τ 1 at (u 1 , λ 1 ) with conserved orientation, i.e., sign τ 0 , τ 1 = 1, can be computed from Although A(s) is singular at (steady) BPs, the convergence of the Newton loop near a branch follows from the Newton-Kantorovich theorem and the so-called Bordering Lemma [17, Lemma 3.1], which deals with the structure of A ∈ R (n+1)×(n+1) , which consists of the main part G u ∈ R n×n and the borders G λ ∈ R n×1 , ξu 0 ∈ R 1×n , and (1 − ξ)λ 0 ∈ R. Similar bordered structures appear again and again in continuation and bifurcation analysis and numerics. Algorithm 2.1 summarizes the basic continuation idea, already including some elementary stepsize ds control.

Algebraic Bifurcation Equations
In the following discussion of numerical branch switching of steady states, (u 0 , λ 0 ) is called a BP, if two or more smooth branches Algorithm 2.1 Branch continuation p=cont(p,aux). Here and in the following the input/output argument p (as in problem) is the pde2path matlab struct which contains all the problem data, and aux stands for auxiliary arguments  (24), set (u 0 , λ 0 , τ 0 )=(u 1 , λ 1 , τ 1 ) and go to 1.

Simple Branch Points
In general, only for m = 1 the QBE determine all (i.e., both) branches through (u 0 , λ 0 ). In this case, (28) reduces to and if (α 0 , α 1 ) is one solution, then the other is distinct (linear independent) if aα 1 + bα 0 = 0. For the branch switching, let (α 0 , α 1 ) with α 0 = λ 0 and α 1 = ψ, u 0 be determined by the branch already computed, and, assuming the generic case Then the other root of (29) is . For normalization we choose α 0 = a, and thus obtain Algorithm 2.2.

Algorithm 2.2
Branch-switching at a simple BP via p=swibra(dir,ptnr), and subsequent continuation of the new branch by p=cont(p). Here the arguments dir,ptnr stand for the pde2path setting that the BP has filename ptnr.mat in directory dir from the branch already computed, and compute a, b, c from (29).
If α 0 = 0, then choose τ = 0 ds as a guess for the bifurcation direction.

Use p=cont(p) to continue the new branch.
Detection of BPs, Exchange of Stability m = 1 in (25) together with aα 1 + bα 0 = 0 are the general bifurcation conditions for simple BPs, also for X a general Banach space, cf. [11]. In our restriction X = R n we further obtain the following results: (B1) [37, §5.8]. If μ 1 (s 0 ) = 0 is an algebraically simple eigenvalue of A, and μ 1 (s) changes sign at s 0 , then (u 0 , λ 0 ) is a simple BP. This yields a simple but efficient criterion to detect BPs:

BDT1 (bifurcation detection test 1).
To detect simple BPs in G : R n+1 → R n we monitor sign changes of det A. This can be done efficiently if we already have an LU-decomposition of A, [37, §5.21]. Under the same assumptions as in (B1), det G u (s) changes sign at s 0 on the branches through (u 0 , λ 0 ) for which μ 1 (s 0 ) = 0.
BDT1 excludes FPs, where a simple eigenvalue of A reaches zero but det A does not change sign. Moreover, BDT1 detects an odd number (counting multiplicities) of eigenvalues crossing, but excludes bifurcations via even numbers of eigenvalues crossing. This excludes Hopf bifurcations, but also, e.g., steady BPs of even multiplicity. Thus, in pde2path we also provide an alternative algorithm BDT2 to detect BPs (and Hopf bifurcation points), which is based on computing a few eigenvalues of G u .
After detection of a bifurcation between s k and s k+1 , the bifurcation can be localized by a bisection method, with a secant, tangent, or quadratic predictor. Although this is a slow method for finding roots of continuous real functions, in the setting of calculating sign changes of det A via LU decomposition (BDT1), and also for BDT2, it seems difficult to improve. Alternatively, see §2.4 for so-called extended systems for fold-, branch-and Hopf-point localization and continuation.
From (B2) we obtain that bifurcations with λ (s 0 ) = 0 from stable branches for Higher Multiplicity The case m ≥ 2 is more difficult, and the QBE (28) may (and typically will) not yield (all) bifurcating branches. The question which order of Taylor expansion in the sense of (27) is needed is called determinacy. In a loose sense, see [87, §6.4] for precise definitions, a given system of algebraic bifurcations equations in α (such as (28)) is called k-determined if any small perturbation of order k + 1 does not qualitatively change the set of (small) solutions. In this sense, transcritical bifurcations are (generically, i.e., unless some special structure occurs) 2-determined, and pitchfork bifurcations are generically 3-determined. Thus, to compute the α for pitchfork branches we need to at least consider the associated cubic bifurcation equations (CBE). In principle, in case of higher order indeterminacies, this must be further continued to higher order, which quickly becomes rather complicated.
In pde2path, we take a practical approach and proceed as follows. The function qswibra searches for solutions (α 0 , α) of the QBE (28) with α 0 = 0 (transcritical bifurcations). 4 All solutions found are stored in p.mat.qtau, and an orthonormal base of the kernel is stored in p.mat.ker, where as in Algorithm 2.2 p is used as the name of the pde2path struct which stores the problem data. Subsequently, we can either select (via seltau) a bifurcating direction from p.mat.qtau, or generate (via gentau) a guess for a bifurcating direction as a linear combination of p.mat.ker, and afterwards call cont.
The approach is summarized in Algorithm 2.3. It is only theoretically sound for 2-determined branches and (in the form cswibra) for 3-determined pitchforks, but it works well and robustly for all the problems we considered so far in pde2path, though some fine tuning (via the optional argument aux, containing tolerances and initial guesses for Newton loops) is sometimes needed. In case of continuous symmetries, further preparatory steps must be taken; see [68] for details.
Unfortunately, there is no general method to detect (31) which can be used for large n u . If (30) comes from a dissipative problem, where most eigenvalues of G u are in the right complex half plane and bounded away from the imaginary axis, then we may try to just compute n eig eigenvalues of G u of smallest modulus, which for moderate n eig can be done efficiently (by inverse vector iteration). Extending this idea, we can also use spectral shifts iω j , j = 1, . . . , n eig , near which we expect eigenvalues to cross the imaginary axis and compute and inspect a few eigenvalues near each iω j . The method, called BDT2 in pde2path, is ad hoc, but with suitable care works quite robustly. If (via BDT2 and, e.g., bisection) we have localized a (simple) HP, then we may want to compute the bifurcating branch of POs, Hopf branch or PO branch in short. Letting μ j (λ) = iω H + O(λ − λ H ), the first order predictor for the bifurcating branch is, from general results on Hopf bifurcations, with step length ds, where ψ is the eigenvector associated to iω H . The continuation of the PO branch is again a predictor-corrector method. First we rescale t = Tt in (18) to obtain with unknown period T , but with initial guess T = 2π/ω H at bifurcation. Since (33) is autonomous (does not explicitly depend on time t), if u H is a PO of (33), then so is any time translate u H (t) = u H (t + δ) of u H . Thus, to obtain unique solutions of (33) we use a PC, for which we choose where ·, · is the scalar product in R n u andu 0 (t) is from the previous continuation step. For the step length condition we choose where again T 0 , λ 0 are from the previous step, ds is the step-length, = d ds denotes differentiation with respect to arclength, ξ H and w T denote weights for the u and T components of the unknown solution, and t 0 = 0 < t 1 < · · · < t m = 1 is the temporal discretization. Thus, the step length is ds in the weighted norm Setting U = (u, T , λ), in each continuation step we need to solve the extended system where G(U ) = 0 is the discretization of (33), including the periodicity condition u m − u 1 = 0. Using Newton's method for (37) we have In (38), ∂ u G ∈ R mn u ×mn u is a large (but sparse) matrix, and A is bordered with borderwidth 2, and bordered elimination solvers such as lssbel [68, §4.5] may yield order of magnitude speedups.

Floquet Multipliers and PO Bifurcations
The (in)stability of a PO u H , and possible bifurcations, are analyzed via the Floquet multipliers γ j , obtained from finding nontrivial solutions (v, γ ) of the variational boundary value problem (in time) The map v(0) → v(1) = Mv(0) defines the so-called monodromy matrix M ∈ R n u ×n u , and the multipliers γ j are the eigenvalues of M. By translational invariance of (33), there always is the trivial multiplier γ 1 = 1, with associated solution v = d dt u H of (39). A simple test for the accuracy of the multiplier computation is the numerical error of the trivial multiplier γ 1 , i.e. err γ 1 := |γ 1 − 1|.
Using the same t-discretization for v in (39) as for u in (33), the multipliers γ j can be computed from the Jacobian A in (38), but this must be done with care, as discussed in, e.g., [28,49], see also [68, §3.4] for more comments, including the pde2path implementation and subsequent PO branch switching.

Further Algorithms and Comments
The algorithms from §2.3 are the basis of arclength continuation of steady states and periodic orbits, including bifurcation detection and localization, and branch switching. Here we comment on two further classes of algorithms, namely extended systems for special point localization and continuation, and PCs, which both play an important role in practical applications of continuation methods.
Extended Systems An alternative to bisection for FP, BP and HP localization is to set up appropriate extended systems. Assume that a FP has been detected between (u, λ)(s a ) and (u, λ)(s b ), i.e., near (u 0 , λ 0 ) = (u, λ)(s 0 ) where for instance we may use just one or two bisection step(s) to approximate (u, λ)(s 0 ) between (u, λ)(s a ) and (u, λ)(s b ). We then set up the extended system so that φ is in the kernel of G u with φ, φ = 1. The bordering lemma [17, Lemma 3.1] yields that the Jacobian is non-singular at a (generic, i.e., quadratic) FP, and hence we may run a Newton loop on (42), using (u 0 , λ 0 ) and the eigenvector φ 0 belonging to the eigenvalue μ 0 of smallest modulus at (u 0 , λ 0 ) as a starting guess, which upon convergence gives us the fold point (u * , λ * ). Similarly, in many problems it is useful to continue FPs in a second parameter to obtain fold curves. Thus freeing a second parameter β, say, and writing w = (λ, β) we set up the extended system where p(U, s) is the analog of the arclength (20), i.e., now

More on Phase Conditions
As noted in Remark 2.7b), in continuation often the PDE M∂ t u = −G(u, λ) must be extended by further equations Q(u, λ) = 0, additional to the arclength condition p(u, λ, s) = 0 as in (19) with p given in (20), modified to ψ(u, λ, T , s) = 0 in (35) for the Hopf case. Two examples of additional equations are already given in (42) for FP continuation, and in (34) as a phase condition (PC) to remove the time-shift of POs from the continuation. Similar PCs are also needed to remove other continuous symmetries, e.g., for spatially translationally invariant problems, which occur for periodic BCs. For instance, over a 1D interval with periodic BCs, a constant coefficient PDE M∂ t u = −G(u) is translationally invariant, and if u with ∂ x u ≡ 0 is a solution of G(u, λ) = 0, so is u(· − ξ) for any ξ ∈ . Hence and ∂ u G always has ∂ x u in its kernel, which is problematic for, e.g., robust convergence of Newton loops, and bifurcation detection. To remove the continuous translational symmetry we add a PC Q(u, λ) = 0, here of the form where u old is either a fixed profile, or the solution from the previous continuation step.
Since ∂ x u old , u old = 1 2´ ∂ x u 2 old dx = 0 we can simplify Q(u, λ) = ∂ x u old , u L 2 , but this makes no difference for the implementation. Moreover, ∂ u Q(u, λ)v = ∂ x u old , v L 2 has a simple form and can be easily provided. The additional equation (47) requires to free an additional parameter η, which may be seen as a Lagrange multiplier for the constraint Q(u, λ) = 0, and we add the generator of the underlying symmetry group to G, i.e., here,

Applications
The MATLAB package pde2path is designed to apply the algorithms from §2, and several more, in a user friendly way to the large class of PDE systems of type (2). Many worked out examples are given in [67,68,76] and in the tutorials at [73], and are included as demos in the software download. See also [16] for a quickstart guide (download and installation instructions), reference card, and an overview of the pde2path demos. 5 Here we apply pde2path to the example problems (5), (6), and (7) and further dead core problems, with the respective demos again available at [73]. These demos follow some standard setup, and here we only give minimal comments on implementations; for details including commented code see the accompanying supplementary information [72].

Turing and Hopf Bifurcations in a 2-Component System
We start with the system (5), i.e., Throughout we fix (τ, d) = (8, 0.05), and initially also α = 0.02, and take j 0 as the primary continuation/bifurcation parameter. For all (j 0 , α), (49) has the unique spatially homogeneous steady state Denoting the terms without derivatives in (49) by f , the linearization of (49) at u * has the form where k ∈ R is called wave number, e ikx φ(k) is called a (Fourier-)mode, and (μ(k), φ(k)) ∈ C × C 2 is an eigenpair ofL(ik) : is called the dispersion relation, where for each k we sort the eigenvalues such that Reμ 2 (k) ≤ Reμ 1 (k) (here we use the more standard sign convention that Reμ > 0 means instability). Over = (−l x , l x ) we additionally have to fulfill the Neumann BCs, which restricts k to the (dual) lattice L = πZ/(2l x ). The solution u * is stable if Reμ 1 (k) < 0 for all k ∈ L, unstable if there exists a k such that Reμ 1 (k) > 0, and instabilities and bifurcations may occur if there is a k = k c such that Reμ 1 (k c ) = 0. The instability is called a long wave instability if k c = 0 and Imμ 1 (k c ) = Imμ 2 (k c ) = 0; it is a Hopf instability if k c = 0 and Imμ 1 (k c ) = Imμ 2 (k c ) = 0, while for k c = 0 (and consequently Imμ 1 (k c ) = 0) it is called a Turing instability. 6 Figure 5(a) shows the Turing and Hopf lines in the j 0 -α plane, where the respective instabilities first occur. These lines can be computed analytically, see [53] (for = R), but here we compute them numerically by BP and HP continuation (see §2.4) on = (−l x , l x ), l x = 8π/k c , with k c = (ατ/D) 1/4 . 7 In Fig. 4(b-d) we show the dispersion relations k → μ 1 (k) at different values of (j 0 , α), where the * illustrate the allowed k over . For given j 0 , the steady state u * is stable for α above the max of the red and blue lines. As we decrease α, (or decreasing j 0 from 3.5 for fixed α = 0.05, say) we either first cross the blue line, meaning that u * looses stability to a Turing mode, or the red line, meaning that u * loses stability to a Hopf mode. Moreover, after crossing into the 'instability region', there are typically further (Hopf or Turing) instabilities in quick succession. There are two codimension-2-points, at In any neighborhood of these, the first instability can be either to a Hopf or to a Turing mode.

Some First Results
To plot BDs for (49), we use the "norm" u 2 defined as (54) for steady states (a), and POs (b), respectively, such that u 2 = 0 means u ≡ u * . This gives better graphical separation of branches in the BDs.
In Fig. 5(a) we show steady state (Turing) branches (T1, blue, to T4, green) bifurcating from the trivial branch u ≡ u * , a secondary steady branch (T1-1 orange), and two secondary Hopf bifurcations (T1-1-h10, red, from the 10th HP on the T1-1 branch, and similarly T1-1-h18 from the 18th HP). 8 The blue branch T1 bifurcates subcritically, and becomes stable in the fold near pt30. The further Turing branches (with sideband wave numbers) behave similarly, and branch T4 goes furthest to the right. Moreover, on all the Turing branches there are secondary bifurcations to localized (steady) patterns, and also some Hopf bifurcations to patterns modulated in time.
Here we first follow the 1st secondary branch (T1-1, orange) which contains (steady) fronts between u = u * and Turing patterns, snaking up by adding one oscillation for every two folds, and with stable segments in direction SE to NW. By even extension of solutions over the left boundary (based on the Neumann BCs), we can also regard such fronts as localized patterns.
First indicated in [57], such snaking branches of localized patterns ("homoclinic snaking", in 1D, 2D and 3D) have seen much interest in recent years, see [2,4,7,10,39,40,74,75], and [68,Ch. 8 and 9] for various examples illustrated with pde2path. In Fig. 5(a), the stable segments of the snake are bounded by folds on the right and HPs on the left, and the red branches are the PO branches bifurcating at the HPs 10 and 18. The solutions on these branches essentially consist of the (steady) Turing pattern on the left, and (small amplitude) oscillations in time around u * on the right. These PO branches connect pairs of HPs on the orange branch, and all only contain unstable solutions, see panel (b2) for two samples. 9 The orange snake connects to the 4th Turing branch T4 near its fold, and thus connects patterns of different wave lengths. Panel (c) gives a partial illustration of the next snake bifurcating from the T1 branch, and one example Hopf branch bifurcating from this second snake, which all again only contain unstable solutions. Finally, there are many further BPs on all the Turing branches, giving a hint of the very rich bifurcation structure of steady states for (49), near the codimension-2 point. Fig. 6 we (re)start with the primary Hopf bifurcation of a spatially homogeneous PO (blue branch) from the trivial branch (for graphical reasons plotted after the 15th point). The branch bifurcates subcritically, and becomes stable at large amplitude after the fold, see, e.g., pt22. The Floquet spectrum at H1/pt16 in (b) illustrates that at lower amplitude there are many unstable directions, and consequently many PO bifurcations before the branch becomes stable. Finally, we consider a mixed mode branch that is not obtained by bifurcation from an already known branch, but by splicing together the Turing pattern T1/pt35 on part of the domain, with the PO from H1/pt25 on the complement. This is just an initial guess for subsequent Newton loops to get on a solution branch, and depending on the choice of, e.g., the "splitpoint" (by trial and error), there may be very large 8 The BDs we show are essentially verbatim outputs of the pde2path scripts (see the plotbra and plotsol commands in cmds1.m in [72]), with minor adjustment of the label placement via click and drag. For instance, the label "30" on the blue branch "T1" means the 30th continuation point. The data is also saved to file accordingly, i.e., as pt30.mat in (sub)directory T1, and we use the same labeling for solution plots (and Floquet multipliers) in (b) and (c). Steady BPs are indicated by •, HPs by , and FPs by ×. 9 One might call these solutions "mixed modes", but following [53] we reserve this term to mixes of Turing patterns and large amplitude Hopf oscillations, which can be stable, see Fig. 6. initial residuals and the Newton loop may fail. However, here we get to pt1 on the red branch, called lh1a, which happens to be stable. 10 Continuing in one direction, the oscillating parts of the solutions shrink, e.g., lh1a/pt40, while in the other direction they grow in a snaking fashion, with stable segments. However, ultimately the branch becomes and stays unstable in both directions, and goes back and forth in j 0 , with solutions loosing more and more structure, and we were not able to determine where it bifurcates from/connects to. This is somewhat unsatisfying, and unfortunately not altogether untypical. Essentially, it means that we have to accept that there are "too many branches" to get a complete picture, but on the other hand Figs. 5 and 6 show that continuation and bifurcation are useful tools to understand big parts of the organization of sets of solutions.

Large Amplitude Mixed Modes In
Comments Where does the multitude of solution branches (of which we only give a small sample in Fig. 5 and Fig. 6) come from? And what do we learn from the continuation/bifurcation approach compared to just finding (stable) solutions (steady or time-periodic) from DNS as in [53]?
Regarding the first question, the dispersion relation(s) Reμ(k) of eigenvalues as a function of wave number k are rather flat. As a consequence, over non-small domains there are many eigenvalues (Turing and/or Hopf) crossing the imaginary axis near criticality in short succession, and there can be sub-harmonic (half wave number or frequency) Turing and/or Hopf modes. This explains the multiplicity of patterns close to bifurcation, which for a similar situation are also analyzed via amplitude equations in [15]; [53] and [15] then use DNS to find a variety of (stable) patterns further away from the primary bifurcation, including some not discussed here, e.g., spiking of subharmonic Turing-Hopf modes. See also [84] for DNS in related systems and a discussion of the role of subcritical bifurcations in localized oscillating patterns. Here we study patterns far from the primary bifurcation via continuation, which for instance allows the computation of the localized steady patterns in the snakes in Fig. 5 in a very efficient way, and which gives further information about how different patterns are connected. Thus, both approaches complement each other. Finally, we may expect a still much richer solution space for (49) in 2D, where already the steady states allow a much richer variety of spots vs stripes (see also §3.2).

Pattern Formation on Disks
Our second example deals with somewhat non-standard patterns due to geometric constraints, and with (spatial) PCs. The SH equation (6), i.e., ∂ t u = −(1 + ) 2 u + εu + νu 3 − u 5 , is a prototypical model for pattern formation. See [63] and, e.g., [12,13,56,61] for reviews. We assume Neumann BCs ∂ n u = ∂ n u = 0 for u and u, and rewrite the 4th order equation (6), as a parabolic-elliptic system, with a singular dynamical mass matrix M d = 1 0 0 0 , f (u 1 ) = νu 3 1 − u 5 1 , and Neumann BCs for u 1 and u 2 . See, e.g., [68,Remark 8.1] for the equivalence of (6) and (55) over convex Lipschitz domains, or general domains with a smooth boundary. We first fix ν > 0, and use ε ∈ R as the primary bifurcation parameter. Over R d , the trivial branch u ≡ 0 loses stability at ε = 0 to modes e i k· x with wave number k = k ≈ 1 and growth rate μ(k) = −(1 − k 2 ) 2 + ε. Essentially the same holds over bounded boxes with Neumann BCs or periodic BCs, with k from the pertinent lattice, depending on the box size. For both, the cubic-quintic case (6) and the quadratic-cubic case with f (u) = νu 2 − u 3 with ν > 0, there are subcritical bifurcations of branches of patterned states, and subsequent secondary bifurcations to snaking branches of localized patterns, and loosely speaking, one difference between the quadratic-cubic and the cubic-quintic case is that the former favors hexagons, while the latter favors stripes. See, e.g., [68,Chap. 8], and the references therein.
Here we report a few results for (55) on a disk, referring to [81] for further details. The SH equation (6), and hence the system (55), is variational with respect to the energy Using integration by parts and the boundary conditions we obtain Hence F decreases along orbits of (6), and since F[u] is bounded from below and since for fixed (ε, ν) we have a discrete set of steady states, every solution converges to a steady state, and in particular there are no Hopf bifurcations and no POs. Introducing polar coordinates (x, y) = ρ(cos φ, sin φ), the primary bifurcation directions (from the trivial branch) are given by Fourier-Bessel modes where m ∈ N 0 is the azimuthal wave number, J m is the mth Bessel function of the first kind, α = 0 or α = π/2, and the constants a m , b m (with a 2 m + b 2 m = 1) and k ± are from a discrete set determined by the boundary conditions. Even for rather small disk radius R, there is a rather large number of instabilities for small ε, see Fig. 7(a), and the order in which they appear depends sensitively on the radius R, see [81].
In a first classification, the modes can be distinguished into axisymmetric ones m = 0 and nonaxisymmetric ones m = 0, and additionally into localized ones (large near the center of the disk) and wall modes (large near the boundary). Following the remarks in §2.4, for the localized modes (axisymmetric or not) it turns out to be useful to switch on translational PCs, while for the non-axisymmetric modes on the full disk we need rotational PCs due to the rotational invariance of the Laplacian and of the disk.
In the following we shall briefly report, for R = 14, on the branches corresponding to the first axisymmetric mode φ 5 and the D − 4 mode φ 6 ( Fig. 7(b)), on some secondary bifurcations from the radial branch (Fig. 8), and on "daisy" branches corresponding to the wall mode φ 3 , and their secondary bifurcations to localized daisies. This focus on the 5th, 6th and 3rd bifurcating branches may appear arbitrary, but these are prototypical branches, and, as already said, the sequence of the bifurcations from u ≡ 0 strongly depends on the radius R. Again, the goal is to obtain a basic understanding of the organization of the set of steady states.

Remark 3.1 a)
Due to the already large number of close-together patterned solutions on rather small disks, potential branch-jumping during continuation becomes a problem, in particular in the neighborhood of bifurcation points. By this we mean the uncontrolled and undetected switching of states from one solution branch to another, typically coming with a loss of symmetry. Obvious ways to mitigate this are choosing a smaller domain (see (b)), and carefully choosing a radially symmetric mesh, see also [68, §4.2]. Additionally, it turns out that choosing a piecewise quadratic (Lagrangian P 2 ) FEM, see [70], instead of the default piecewise linear P 1 FEM helps to mitigate branch jumping.  [81] and the tutorial at [80] we present our results only on the half disk, where also at x 1 = 0 we impose Neumann BCs for (u 1 , u 2 ). The solutions can then be mirrored to the full disk. This simplifies the numerics, as we need less degrees of freedom for the discretization, and as for localized solutions (see below) we do not need a PC in x 1 , and, equally importantly, keeping symmetry of solution branches is easier on the smaller domain. However, the half disk also reduces the number of allowed modes to α = 0 in (58). Therefore, when below discussing stability of nontrivial solutions, we always first extend them to the full disk and then compute the spectrum.
Some Results In Fig. 7(b) we present a bifurcation diagram of the branches from φ 5 and φ 6 , for ν = 2. 11 The black branch bifurcates at ε = ε 5 ≈ 0.88 × 10 −2 in direction φ 5 , and contains axisymmetric states which start as a spot at small norm; the branch then snakes in ε, containing short stable segments, before turning into a domain filling stable target state in the last fold. However, this behavior depends in a very sensitive and interesting way on the disk radius. For instance, for R = 15 the primary spot branch does not connect to a target state, but instead reconnects to u = 0 at slightly larger ε, as discussed in detail in [81]. The red branch in Fig. 7(b) is an example of states with D − m symmetry, which means that the states are invariant under rotation by 2π/m together with a sign change u → −u. The branch starts with localized D − 4 states (hence again needs a y-PC); along the branch, the arms first extend towards the wall (state 90), and then broaden, and the SW → NE segments of this branch contain stable segments.  Figure 8 shows some secondary bifurcations from the radial branch, starting with 2-arms, 3-arms and 4-arms branches (with D m symmetry, i.e., invariance under rotation by 2π/m) bifurcating at small norm (a-c). Again the arms grow towards the wall, and the 3-arms and 4-arms branches broaden as they reach the wall, and seem to reconnect to the axisymmetric black branch near its upper left fold. The 2-arms branch in (a) behaves differently and rather loses spots after point 4. In (d) we look at the reconnection more closely. It turns out that the 3-arms and 4-arms branches do not connect directly to the target, but to a "crown" state which bifurcates from the radial branch near its last fold, and has D 12 angular structure at the wall; see sample 5 in Fig. 8(d). Conversely, there are 1-arm and 2-arm states bifurcating from the crown branch; see the critical eigenvectors at BPs 1 and 2 in (d), and the further discussion of these bifurcating branches in [81].
In Fig. 9(a) we set ν = 1.4 and consider the blue branch bifurcating at ε = ε 3 ≈ 0.292 × 10 −3 in direction of the wall mode φ 3 . Following [46] we call these solutions daisy states. At the first BP on the daisy branch, a branch (green) of localized daisies bifurcates, snakes up by adding a petal at every other fold, and reconnects to the blue daisy branch near its top left fold, where the daisies become stable. This is very similar to the snaking of the T1-1 branch in Fig. 5, or to the snaking of localized periodic patterns in the 1D SH equation [8]; like there, here we also have an intertwined 2nd branch of odd localized daisies, and the two are connected by rungs, cf. [81, Fig. 12].
However, the quasi 1D snaking on the R = 14 disk only holds at moderate subcriticality, i.e., moderate ν. For instance, for ν = 2 there also bifurcates a localized 1-petal daisy branch from the daisy branch at low norm, and a 1-hole "plucked" daisy branch near the fold, but the two do not connect. To understand this breakup, in [81] we use some fold continuation in ν of the folds as in Fig. 9(a). It turns out that these folds do not continue to ν = 2, but the fold continuation folds back at some ν < 2. In Fig. 9(b) we continued FP12 from (a) to ν = 1.5, and then switch back to continuation in ε. In the "up" direction, the branch (in green) behaves as in (a), but in "down" direction the branch (in red) does not loose further petals to connect to the daisy branch at low norm. Instead, the petals start to grow towards the center of the domain. Such states of patches of stripes are also known for the SH equation in the plane (without boundary) and called worms. Thus, in [81] we term states like 100 "boundary worms".
Here we end the review of [81], which discusses the above results and a number of further patterns in more detail. One crucial point is that continuation and bifurcation software helps to uncover the very rich solution structure in a more systematic way than just DNS. All of the above branches (and many more in [81]) start out unstable, but later contain stable segments, and (secondary, tertiary,..) bifurcations to further branches with stable parts.

Experiments on Dead Core Pattern Formation
Our last (class of) example(s) is "experimental", in the sense that the standard analytical bifurcation theory does not apply because the problems are not smooth, and not even piecewise Lipschitz. This also reflects in the need to modify some of the pde2path algorithms.
A classical example of a dead core (DC) problem is given by the elliptic equation where ⊂ R d is a bounded domain, λ > 0 a parameter, and f : R → R a function such that The associated parabolic problem is Setting v = u γ , (61) with f (u) = u γ + , u + = max(u, 0), is equivalent to the degenerate diffusion problem where f (v δ ) = v, which however does not seem simpler to analyze than (61). For large λ, (59) has DC solutions which feature a subdomain 0 ⊂ where u ≡ 0. The motivation comes from reaction-diffusion problems where the reaction is so fast near u = 0 that the diffusion cannot supply further reactants and a DC without reactants develops. For (59) in 1D with f (u) = u γ + and wlog = (0, 1), DC solutions can be computed explicitly. For instance, for γ = 1/2 the ansatz yields β = 4, α = λ 2 /144 and δ = 1/2 − √ 12/λ, i.e., a DC for λ > 48, which turns out to be the unique non-negative solution for λ > 48, while for λ < 48 the unique solution is positive. The DC grows as λ → ∞, and the solution is a classical one and in fact in C 3 ([0, 1]), while the positive solutions for λ < 48 are C ∞ . Similar computations can also be done for balls in d = 2, 3 (2D, 3D cases), and are useful for numerical checks.
Problems of this kind, in d space dimensions and under various BCs and general assumptions, have been studied analytically in detail, including also the case of degenerate diffusion and quasilinear cases. We refer to [29] and the references therein for various older results on (59), which can be summarized as follows: DC solutions exist for many problem classes, and are generally rather smooth, i.e., at least classical solutions; the dead cores 0 = {x ∈ : u(x) = 0} are convex if is convex, and their "free boundaries" ∂ 0 determined together with the solution itself are C 2+γ if ∂ is smooth. See also [22,47,64] for more recent results and, e.g., [22] for a comprehensive review of the physical background including the case of systems of equations, for which we also refer to, e.g., [18].
For (59), some results on multiplicity of solutions are available, see for instance [30] and the references therein. These yield that for small λ (relative to the domain size) there is a unique positive solution, and for λ sufficiently large there is a unique DC solution, while [30] also discusses uniqueness vs non-uniqueness in DC problems with BCs depending on a parameter λ. However, a genuine bifurcation analysis of DC solutions seems to be largely lacking. The existing theory for bifurcations in nonsmooth dynamical systems, e.g., [20,21,48,50], does not apply as it deals with at least piecewise Lipschitz systems, with the notable exception of [43] where scaling laws near folds in non-Lipschitz systems are derived.

Numerical Issues and Implementation
In the demo acdc, available at [73], we numerically study slight modifications of (59) in 1D and 2D, which (via x-dependent nonlinearities [1]) allow various stable branches with multiple interfaces bifurcating from u ≡ 1, which subsequently develop different DCs. Here we give a few general comments on modifications of pde2path needed for DC problems, essentially sticking with (59), which then can immediately be transferred to the DC Schnakenberg system studied below.
The convergence theory for the P 1 -FEM for equations of type (59) is discussed in, e.g., [54] and [3]. The basic idea of [54] is to apply the FEM to regularized problems where f is replaced by globally Lipschitz regularizations f ε where ε = Ch 2 depends on the FEM discretization parameter h (mesh-width). In a nutshell, this yields the a-priori estimate u − u ε,h ∞ ≤ C(h log(1/h)) 2 for h < h 0 . In [3], this is combined with a nonlinear successive-over-relaxation method to solve the (nonsmooth, at u = 0) nonlinear algebraic systems obtained from the FEM applied to the original equation (59). This uses either trivial (linear) steps, or Newton steps, or regula falsi steps, depending on the residual, and the method is shown to be globally convergent to the (unique, due to monotonicity of f (u) = u γ + ) solution of the discrete problem.
Here we aim to use the standard setup of pde2path, with small modifications of the Newton loops u n+1 = u n − (∂ u G(u n )) −1 G(u n ) to ensure that u ≥ 0 everywhere, and we need to deal with the non-differentiability of f at u = 0. To "approximate" ∂ u G(u) (even if it does not exist) in a simple way, we can use forward finite differences, i.e., numerical Jacobians. Alternatively, we implement an approximate ∂ u G by modifying ∂ u (u γ + ) = γ u γ −1 + to γ max(u + , δ) γ −1 with a small δ (e.g., δ = 10 −12 ). Thus, instead of ∞ we have some large number in ∂ u G where u = 0. Adapting the Newton loops as usual works by modifying local (in the problem directory) copies of the respective library functions. For instance, in a local copy of nloop.m we add the command u1(1:p.nu)=max(u1(1:p.nu),0) after computing the update u1 = u n+1 , see [72] for details. The "approximate" ∂ u G together with the modified Newton loops then yields rather robust and fast convergence even in the presence of a DC, and since we do not modify G itself we solve the original singular problem. 12 Of course, if a DC has formed, we then also have a wrong (numerical) ∂ u G for stability and bifurcation computations, and the consequences of this still need to be studied. From the experimental point of view, we find the following hypotheses: H1. The stability information obtained from the eigenvalues of the "approximate" ∂ u G as formed above appears correct in the sense that it fully agrees with DNS. Moreover, if there are unstable eigenvalues, the associated eigenvectors have supports outside the DC. H2. If bifurcations from states with a DC occur (notably in the Schnakenberg problem studied next), then the bifurcation directions have supports outside the DC.
Dead Core Schnakenberg The (classical) Schnakenberg model is a 2-species predator-prey model for cubic autocatalysis: substance u decays at rate 1, but catalyzes its own production while preying on v, which is fed at some constant rate λ. The diffusion constants of u and v are 1 and d 1, respectively, and in summary a basic model is usually considered with homogeneous Neumann BCs. The homogeneous steady state (u, v)=(λ, 1/λ) is stable for large λ, but for d 1 shows a Turing instability for λ decreasing below some (explicitly computable) λ 0 . The version (64) has been a model problem for pde2path since [76]. See also [74] and [68,Ch.9] for numerical bifurcation analysis for (64), in 1D, 2D and 3D, including branches of localized patterns within patterns, for instance snaking branches of localized hexagon patches embedded in stripes, and, e.g., [41,66] for further recent results.
Here we modify (64) to such that u 1/m + becomes singular at u=0 for m>1. The homogeneous branch now is (u, v)=(λ m , λ 1−2m ). Its Turing (in)stability analysis is completely analogous to (64) (m = 1). Again we find Turing bifurcations for decreasing λ (d = 100, say, and then λ in the range 1 to 3, say, such that u hom = λ m is far from the singular point u = 0). Here we are mostly interested in the DCs formed by Turing patterns, and Hopf bifurcations from Turing branches with DCs at low λ.
Numerically, we follow the procedures explained above for (61), see [72] for details. Figure 10 shows results for (65) in 1D. The BD in (a) shows the first three Turing branches T1,T2 and T3 (blue and green), with 2, 2.5, and 1.5 waves in the domain, respectively. These are initially quite similar to the regular case, but T1 and T3, which have longer wavelengths than T2, at low λ develop DCs, which we numerically define as u < 10 −8 and indicate by magenta crosses in (b). Additionally we show two secondary branches T1-1 (red, connecting the first BPs on T1 and T3) and T3-2 (brown), bifurcating at the 2nd BP from T3, where there already is a DC on T3. The critical eigenvector is localized at the right bump of T3/bpt2, and the bifurcating branch decreases this bump.
At yet lower λ we find HPs on T1, T2 and T3, and on T3 we have two DCs at the HP, while on T1 and T3 we still have u 1 > 10 −8 in the valleys. In all cases, the Hopf eigenvectors are localized at the bumps, and consequently in the bifurcating POs mostly the bumps oscillate. See (c), where the time series are from the right boundary point. For T3-H1, we initially keep the DCs and the oscillations are only in the "live" part u > 0, and even at large amplitude in T3-H1/pt8, |u(12, t)| < 10 −12 most of the time. Figure 11(b) shows two DC Turing-patterns in 2D. The domain is = (−l x , l x ) × (−l x / √ 3, l x / √ 3) which yields a hexagonal dual lattice [68, §2.5], and hence double BPs at λ c . In pde2path, such BPs of higher multiplicity can be handled by numerically deriving and solving the quadratic bifurcation equations (28), see Algorithm 2.2.
Here we have two branches at the primary bifurcation, namely stripes (1D solutions extended homogeneously in y-direction, bifurcating in a pitchfork), and (hexagonal) spots. In Fig. 11 we focus on the (transcritical) hexagon branch, and find that in both directions it develops DCs as essentially predicted at bifurcation. As long as u > 0 in (65) the system altogether behaves similar to the regular case (m = 1), and hence we can also find various sorts of patterns as described in [74], which as λ decreases yield various sorts of DCs.
Funding Note Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.