# Professor Gianne Derks

## Academic and research departments

Department of Mathematics, Centre for Mathematical and Computational Biology.### Biography

I graduated as an "Ingenieur in Mathematics" at the Technical University of Eindhoven in 1988 and I got my PhD from the University of Twente in 1992. In 1993/94 I was awarded a NATO Research Fellowship to do research at the University of California Santa Cruz and the MSRI in Berkeley. In the following year 1994/95, I was a Research Fellow at the Simon Fraser University in Vancouver. In 1995, I started at the University of Surrey and I am currently a Professor in the Department of Mathematics.

### Research

### Research interests

- Hamiltonian systems (ODEs and PDEs) with perturbations like dissipation and/or forcing;
- Dissipation induced instability;
- Inhomogeneous wave equations;
- Stability and instability of travelling waves and front solutions, for example for (semi)-kinks in Josephson junctions or ionization waves in gas discharges;
- Multi-symplectic systems;
- Mathematical models in pharmacology, toxicology and metabolism (including pharmacokinetics-pharmacodynamics and physiology based dynamical models);
- Cycles in sleep-wake models, endocrinology;
- Wave dynamics in DNA-RNAP interactions;
- I'm a member of the LMS sponsored Network on Applied Geometric Mechanics;
- I'm a member of the LMS sponsored Network on Mathematics in Life Sciences;
- I'm on the organising committee of the Quantitative Systems Pharmacology-UK (QSP-UK) network;
- I'm on the editorial board of the Journal of Geometric Mechanics.

Further details can be found on my personal web page.

A full list of my publications and pre-prints can be found on my personal web page.

### My teaching

In 20/21 I am lecturing the following module:

- MAT3046 Game Theory for level 3 Mathematical Studies students in Semester 2 (jointly with Prof Anne Skeldon).

Take a look at some suggestions for final year projects or literature reviews (PDF).

### My publications

### Publications

Sleep is essential for the maintenance of human life, yet many features of sleep are poorly understood and mathematical models are an important tool for probing proposed biological mechanisms. The most well-known mathematical model of sleep regulation, the two-process model, models the sleep-wake cycle by two oscillators: a circadian oscillator and a homeostatic oscillator. An alternative, more recent, model considers the reciprocal interaction of sleep promoting neurons and the ascending arousal system regulated by homeostatic and circadian processes. Here we show there are fundamental similarities between these two models. The implications are illustrated with two important sleep-wake phenomena. Firstly, we show that in the two-process model, transitions between different numbers of daily sleep episodes can be classified as grazing bifurcations. This provides the theoretical underpinning for numerical results showing that the sleep patterns of many mammals can be explained by the reciprocal interaction model. Secondly, we show that when sleep deprivation disrupts the sleep-wake cycle, ostensibly different measures of sleepiness in the two models are closely related. The demonstration of the mathematical similarities of the two models is important because not only does it it allow some features of the two-process model to be interpreted physiologically but it also means that knowledge gained from the study of the two-process model can be used to inform understanding of the behaviour of the mutual inhibition model. This is important because the mutual inhibition model and its extensions are increasingly being used as a tool to understand a diverse range of sleep-wake phenomena sucah as the design of optimal shift-patterns, yet the values it uses for the parameters associated with the circadian and homeostatic processes are very different from those that have been experimentally measured in the context of the two-process model

Nitisinone or 2-(2-nitro-4-trifluoromethylbenzoyl)cyclohexane-1,3-dione, is a reversible inhibitor of 4- hydroxyphenylpyruvate dioxygenase (HPPD), an enzyme important in tyrosine catabolism. Today, nitisinone is successfully used to treat Hereditary Tyrosinaemia type 1, although its original expected role was as a herbicide. In laboratory animals, treatment with nitisinone leads to the elevation of plasma tyrosine (tyrosinaemia). In rats and Beagle dogs, repeat low-dose exposure to nitisinone leads to corneal opacities whilst similar studies in the mouse and Rhesus monkey showed no comparable toxicities or other treatment related findings. The differences in toxicological sensitivities have been related to the upper limit of the concentration of tyrosine that accumulates in plasma, which is driven by the amount/activity of tyrosine aminotransferase. A physiologically based, pharmacodynamics ordinary differential equation model of HPPD inhibition to bolus exposure of nitisinone in vivo is presented. Going beyond traditional approaches, asymptotic analysis is used to separate the different timescales of events involved in HPPD inhibition and tyrosinaemia. This analysis elucidates, in terms of the model parameters, a critical inhibitor concentration (at which tyrosine concentration starts to rise) and highlights the contribution of in vitro measured parameters to events in an in vivo system. Furthermore, using parameter-fitting methods, a systematically derived reduced model is shown to fit well to rat data, making explicit how the parameters are informed by such data. This model in combination with in vitro descriptors has potential as a surrogate for animal experimentation to predict tyrosinaemia, and further development can extend its application to other related medical scenarios.

We consider inhomogeneous non-linear wave equations of the type u =u +V (u, x)-αu (α≥0). The spatial real axis is divided in intervals I , i=0,..., N+1 and on each individual interval the potential is homogeneous, i.e., V(u, x)=V (u) for x∈I . By varying the lengths of the middle intervals, typically one can obtain large families of stationary front or solitary wave solutions. In these families, the lengths are functions of the energies associated with the potentials V . In this paper we show that the existence of an eigenvalue zero of the linearisation operator about such a front or stationary wave is related to zeroes of the determinant of a Jacobian associated to the length functions. Furthermore, the methods by which the result is obtained is fully constructive and can subsequently be used to deduce the stability and instability of stationary fronts or solitary waves, as will be illustrated in examples. © 2012 Elsevier Inc.

Operators on unbounded domains may acquire eigenvalues that are embedded in the essential spectrum. Determining the fate of these embedded eigenvalues under small perturbations of the underlying operator is a challenging task, and the persistence properties of such eigenvalues are linked intimately to the multiplicity of the essential spectrum. In this paper, we consider the planar bilaplacian with potential and show that the set of potentials for which an embedded eigenvalue persists is locally an infinite-dimensional manifold with infinite codimension in an appropriate space of potentials.

A long Josephson junction containing regions with a phase shift of is considered. By exploiting the defect modes due to the discontinuities present in the system, it is shown that Josephson junctions with phase shift can be an ideal setting for studying localized mode interactions. A phase-shift configuration acting as a double-well potential is considered and shown to admit mode tunnelings between the wells. When the phase-shift configuration is periodic, it is shown that localized excitations forming bright and dark solitons can be created. Multimode approximations are derived confirming the numerical results.

*SIAM Journal on Mathematical Analysis*33pp. 1356-1378

The linear stability problem for solitary wave states of the Kawahara---or fifth-order KdV-type---equation and its generalizations is considered. A new formulation of the stability problem in terms of the symplectic Evans matrix is presented. The formulation is based on a multisymplectification of the Kawahara equation, and leads to a new characterization of the basic solitary wave, including changes in the state at infinity represented by embedding the solitary wave in a multiparameter family. The theory is used to give a rigorous geometric sufficient condition for instability. The theory is abstract and applies to a wide range of solitary wave states. For example, the theory is applied to the families of solitary waves found by Kichenassamy--Olver and Levandosky.

*φ*

_{0}/2-flux case, In: Physical Review B69(212503)

We consider Josephson vortices at tricrystal boundaries. We discuss the specific case of a tricrystal boundary with a pi junction as one of the three arms. It is recently shown that the static system admits an (n+1/2) phi0 flux, n=0,1,2 [Phys. Rev. B 61, 9122 (2000)]. Here we present an analysis to calculate the linear stability of the admitted states. In particular, we calculate the stability of a 3 phi0/2 flux. This state is of interest, since energetically this state is preferable for some combinations of Josephson lengths, but we show that in general it is linearly unstable. Finally, we propose a system that can have a stable (n+1/2)phi0 state.

The linear stability of solitary-wave or front solutions of Hamiltonian evolutionary equations, which are equivariant with respect to a Lie group, is studied. The organizing centre for the analysis is a multisymplectic formulation of Hamiltonian PDEs where a distinct symplectic operator is assigned for time and space. This separation of symplectic structures leads to new characterizations of the following components of the analysis. The states at infinity are characterized as manifolds of relative equilibria associated with the spatial symplectic structure. The momentum of the connecting orbit, or shape of the solitary wave, considered as a heteroclinic orbit in a phase-space representation, is given a new characterization as a one-form on the tangent space to the heteroclinic manifold and this one-form is a restriction of the temporal symplectic structure. For the linear stability analysis, a new symplectic characterization of the Evans function and its derivatives are obtained, leading to an abstract geometric proof of instability for a large class of solitary-wave states of equivariant Hamiltonian evolutionary PDEs. The theory sheds new light on several well-known models: the gKdV equation, a Boussinesq system and a nonlinear wave equation. The generalization to solitary waves associated with multidimensional heteroclinic manifolds and the implications for solitary waves or fronts which are biasymptotic to invariant manifolds such as periodic states are also discussed.

Motivated by a problem from pharmacology, we consider a general two parameter slow-fast system in which the critical set consists of a one dimensional manifold and a two dimensional manifold, intersecting transversally at the origin. Using geometric desingularisation, we show that for a subset of the parameter set there is an exchange of stabilities between the attracting components of the critical set and the direction of the continuation can be expressed in terms of the parameters.

Hamiltonian evolution equations which are equivariant with respect to the action of a Lie group are models for physical phenomena such as oceanographic flows, optical fibres and atmospheric flows, and such systems often have a wide variety of solitary wave or front solutions. In this paper, we present a new symplectic framework for analyzing the spectral problem associated with the linearization about such solitary waves and fronts. At the heart of the analysis is a multi-symplectic formulation of Hamiltonian partial differential equations where a distinct symplectic structure is assigned for the time and space directions, with a third symplectic structure -- with two-form denoted by Omega - associated with a coordinate frame moving at the speed of the wave. This leads to a geometric decomposition and symplectification of the Evans function formulation for the linearization about solitary waves and fronts. We introduce the concept of the "symplectic Evans matrix", a matrix consisting of restricted "Omega-symplectic" forms. By applying Hodge duality to the exterior algebra formulation of the Evans function, we find that the zeros of the Evans function correspond to zeros of the determinant of the symplectic Evans matrix. Based on this formulation, we prove several new properties of the Evans function. Restricting the spectral parameter lambda to the real axis, we obtain rigorous results on the derivatives of the Evans function near the origin, based solely on the abstract geometry of the equations, and results for the large $|____lambda|$ behaviour which use primarily the symplectic structure, but also extend to the non-symplectic case. The Lie group symmetry affects the Evans function by generating zero eigenvalues of large multiplicity in the so-called systems at infinity. We present a new geometric theory which describes precisely how these zero eigenvalues behave under perturbation. By combining all these results, a new rigorous sufficient condition for instability of solitary waves and fronts is obtained. The theory applies to a large class of solitary waves and fronts including waves which are biasymptotic to a nonconstant manifold of states as $|x|$ tends to infinity. To illustrate the theory, it is applied to three examples: a Boussinesq model from oceanography, a class of nonlinear Schrodinger equations from optics and a nonlinear Klein-Gordon equation from atmospheric dynamics.

The spectral problem associated with the linearization about solitary waves of spinor systems or optical coupled mode equations supporting gap solitons is formulated in terms of the Evans function, a complex analytic function whose zeros correspond to eigenvalues. These problems may exhibit oscillatory instabilities where eigenvalues detach from the edges of the continuous spectrum - so-called edge bifurcations. A numerical framework, based on a fast robust shooting algorithm using exterior algebra, is described. The complete algorithm is robust in the sense that it does not produce spurious unstable eigenvalues. The algorithm allows us to locate exactly where the unstable discrete eigenvalues detach from the continuous spectrum. Moreover, the algorithm allows for stable shooting along multidimensional stable and unstable manifolds. The method is illustrated by computing the stability and instability of gap solitary waves of a coupled mode model.

We consider the possibility of free receptor (antigen/cytokine) levels rebounding to higher than the baseline level after the application of an antibody drug using a target-mediated drug disposition model. It is assumed that the receptor synthesis rate experiences homeostatic feedback from the receptor levels. It is shown for a very fast feedback response, that the occurrence of rebound is determined by the ratio of the elimination rates, in a very similar way as for no feedback. However, for a slow feedback response, there will always be rebound. This result is illustrated with an example involving the drug efalizumab for patients with psoriasis. It is shown that slow feedback can be a plausible explanation for the observed rebound in this example.

We consider Josephson vortices at tricrystal boundaries. We discuss the specific case of a tricrystal boundary with a pi junction as one of the three arms. It is recently shown that the static system admits an (n+1/2)phi(0) flux, n=0,1,2 [ Phys. Rev. B 61, 9122 (2000) ]. Here we present an analysis to calculate the linear stability of the admitted states. In particular, we calculate the stability of a 3phi(0)/2 flux. This state is of interest, since energetically this state is preferable for some combinations of Josephson lengths, but we show that in general it is linearly unstable. Finally, we propose a system that can have a stable (n+1/2)phi(0) state.

We investigate the existence of stationary fronts in a coupled system of two sine-Gordon equations with a smooth, “hat-like” spatial inhomogeneity. The spatial inhomogeneity corresponds to a spatially dependent scaling of the sine-Gordon potential term. The uncoupled inhomogeneous sine-Gordon equation has stable stationary front solutions that persist in the coupled system. Carrying out a numerical investigation it is found that these inhomogeneous sine-Gordon fronts loose stability, provided the coupling between the two inhomogeneous sine-Gordon equations is strong enough, with new stable fronts bifurcating. In order to analytically study the bifurcating fronts, we first approximate the smooth spatial inhomogeneity by a piecewise constant function. With this approximation, we prove analytically the existence of a pitchfork bifurcation. To complete the argument, we prove that transverse fronts for a piecewise constant inhomogeneity persist for the smooth “hat-like” spatial inhomogeneity by introducing a fast-slow structure and using geometric singular perturbation theory.

Families of stable stationary solutions of the two-dimensional incompressible homogeneous Euler and ideal reduced magnetohydrodynamic equations are shown to be attracting for the corresponding viscous perturbations of these systems, i.e. for the Navier-Stokes and the reduced viscous MHD equations with magnetic diffusion. Each solution curve of the dissipative system starting in a cone around the family of stationary solutions of the unperturbed conservative system defines a shadowing curve which attracts the dissipative solution in an exponential manner. As a consequence, the specific exponential decay rates are also determined. The techniques to analyse these two equations can be applied to other dissipative perturbations of Hamiltonian systems. The method in its general setting is also presented.

We consider a spatially nonautonomous discrete sine-Gordon equation with constant forcing and its continuum limit(s) to model a 0-pi Josephson junction with an applied bias current. The continuum limits correspond to the strong coupling limit of the discrete system. The nonautonomous character is due to the presence of a discontinuity point, namely, a jump of pi in the sine- Gordon phase. The continuum model admits static solitary waves which are called pi-kinks and are attached to the discontinuity point. For small forcing, there are three types of pi-kinks. We show that one of the kinks is stable and the others are unstable. There is a critical value of the forcing beyond which all static pi-kinks fail to exist. Up to this value, the (in)stability of the pi-kinks can be established analytically in the strong coupling limits. Applying a forcing above the critical value causes the nucleation of 2pi-kinks and -antikinks. Besides a pi-kink, the unforced system also admits a static 3pi-kink. This state is unstable in the continuum models. By combining analytical and numerical methods in the discrete model, it is shown that the stable pi-kink remains stable and that the unstable pi-kinks cannot be stabilized by decreasing the coupling. The 3pi-kink does become stable in the discrete model when the coupling is sufficiently weak.

We consider a spatially nonautonomous discrete sine-Gordon equation with constant forcing and its continuum limit(s) to model a 0-pi Josephson junction with an applied bias current. The continuum limits correspond to the strong coupling limit of the discrete system. The nonautonomous character is due to the presence of a discontinuity point, namely, a jump of pi in the sine-Gordon phase. The continuum model admits static solitary waves which are called pi-kinks and are attached to the discontinuity point. For small forcing, there are three types of pi-kinks. We show that one of the kinks is stable and the others are unstable. There is a critical value of the forcing beyond which all static pi- kinks fail to exist. Up to this value, the (in) stability of the pi-kinks can be established analytically in the strong coupling limits. Applying a forcing above the critical value causes the nucleation of 2 pi-kinks and -antikinks. Besides a pi- kink, the unforced system also admits a static 3 pi-kink. This state is unstable in the continuum models. By combining analytical and numerical methods in the discrete model, it is shown that the stable pi-kink remains stable and that the unstable pi-kinks cannot be stabilized by decreasing the coupling. The 3 pi- kink does become stable in the discrete model when the coupling is sufficiently weak.

In this paper the non-linear wave equation with a spatial inhomogeneity is considered. The inhomogeneity splits the unbounded spatial domain into three or more intervals, on each of which the non-linear wave equation is homogeneous. In such setting, there often exist multiple stationary fronts. In this paper we present a necessary and sufficient stability criterion in terms of the length of the middle interval(s) and the energy associated with the front in these interval(s). To prove this criterion, it is shown that critical points of the length function and zeros of the linearisation have the same order. Furthermore, the Evans function is used to identify the stable branch. The criterion is illustrated with an example which shows the existence of bi-stability: two stable fronts, one of which is non-monotonic. The Evans function also give a sufficient instability criterion in terms of the derivative of the length function.

This paper presents an introduction to the existence and stability of stationary fronts in wave equations with finite length spatial inhomogeneities. The main focus will be on wave equations with one or two inhomogeneities. It will be shown that the fronts come in families. The front solutions provide a parameterisation of the length of the inhomogeneities in terms of the local energy of the potential in the inhomogeneity. The stability of the fronts is determined by analysing (constrained) critical points of those length functions. Amongst others, it will shown that inhomogeneities can stabilise non-monotonic fronts. Furthermore it is demonstrated that bi-stability can occur in such systems.

Sleep is essential for the maintenance of human life, yet many features of sleep are poorly understood and mathematical models are an important tool for probing proposed biological mechanisms. The most well-known mathematical model of sleep regulation, the two-process model, models the sleep-wake cycle by two oscillators: a circadian oscillator and a homeostatic oscillator. An alternative, more recent, model considers the reciprocal interaction of sleep promoting neurons and the ascending arousal system regulated by homeostatic and circadian processes. Here we show there are fundamental similarities between these two models. The implications are illustrated with two important sleep-wake phenomena. Firstly, we show that in the two-process model, transitions between different numbers of daily sleep episodes can be classified as grazing bifurcations. This provides the theoretical underpinning for numerical results showing that the sleep patterns of many mammals can be explained by the reciprocal interaction model. Secondly, we show that when sleep deprivation disrupts the sleep-wake cycle, ostensibly different measures of sleepiness in the two models are closely related. The demonstration of the mathematical similarities of the two models is important because it brings a deeper insight to sleep-wake dynamics and allows some features of the two-process model to be interpreted physiologically.

In part I, the Two-Process model for sleep-wake regulation was discussed and it was shown that it could usefully be represented as a one-dimensional map with discontinuities. Here we discuss some recent, more physiological, models of sleep wake dynamics. We describe how their fast-slow structure means that one can expect them to inherit many of the dynamical features of the Two-Process model.

For more than thirty years the `two process model' has played a central role in the under- standing of sleep/wake regulation. This ostensibly simple model is an interesting example of a nonsmooth dynamical system whose rich dynamical structure has been relatively un- explored. The two process model can be framed as a one-dimensional map of the circle which, for some parameter regimes, has gaps. We show how border collision bifurcations that arise naturally in maps with gaps extend and supplement the Arnold tongue saddle- node bifurcation set that is a feature of continuous circle maps. The novel picture that results shows how the periodic solutions that are created by saddle-node bifurcations in continuous maps transition to periodic solutions created by period-adding bifurcations as seen in maps with gaps.

Sleep-wake regulation is an example of a system with multiple timescales, with switching between sleep and wake states occurring in minutes but the states of wake or sleep usually existing for some hours. Here we discuss some general features of models of sleep-wake regulation. We show that some typical models of sleep-wake regulation can be reduced to one-dimensional maps with discontinuities, and show that this reduction is useful in understanding some of the dynamical behaviour seen in sleep-wake models.

The uniformly damped Korteweg--de Vries (KdV) equation with periodic boundary conditions can be viewed as a Hamiltonian system with dissipation added. The KdV equation is the Hamiltonian part and it has a two-dimensional family of relative equilibria. These relative equilibria are space-periodic soliton-like waves, known as cnoidal waves. Solutions of the dissipative system, starting near a cnoidal wave, are approximated with a long curve on the family of cnoidal waves. This approximation curve consists of a quasi-static succession of cnoidal waves. The approximation process is sharp in the sense that as a solution tends to zero as t to infinity, the difference between the solution and the approximation tends to zero in a norm that sharply picks out their difference in shape. More explicitly, the difference in shape between a solution and a quasi-static cnoidal-wave approximation is of the order of the damping rate times the norm of the cnoidal-wave at each instant.

Streamer ionization fronts are pulled fronts that propagate into a linearly unstable state; the spatial decay of the initial condition of a planar front selects dynamically one specific long-time attractor out of a continuous family. A stability analysis for perturbations in the transverse direction has to take these features into account. In this paper we show how to apply the Evans function in a weighted space for this stability analysis. Zeros of the Evans function indicate the intersection of the stable and unstable manifolds; they are used to determine the eigenvalues. Within this Evans function framework, we define a numerical dynamical systems method for the calculation of the dispersion relation as an eigenvalue problem. We also derive dispersion curves for different values of the electron diffusion constant and of the electric field ahead of the front. Numerical solutions of the initial value problem confirm the eigenvalue calculations. The numerical work is complemented with an analysis of the Evans function leading to analytical expressions for the dispersion relation in the limit of small and large wave numbers. The paper concludes with a fit formula for intermediate wave numbers. This empirical fit supports the conjecture that the smallest unstable wave length of the Laplacian instability is proportional to the diffusion length that characterizes the leading edge of the pulled ionization front. © 2008 Springer Science+Business Media, LLC.

We consider a Josephson junction system installed with a finite length inhomogeneity, either of micro-resistor or micro-resonator type. The system can be modelled by a sine-Gordon equation with a piecewise-constant function to represent the varying Josephson tunneling critical current. The existence of pinned fluxons depends on the length of the inhomogeneity, the variation in the Josephson tunneling critical current and the applied bias current. We establish that a system may either not be able to sustain a pinned fluxon, or - for instance by varying the length of the inhomogeneity - may exhibit various different types of pinned fluxons. Our stability analysis shows that changes of stability can only occur at critical points of the length of the inhomogeneity as a function of the (Hamiltonian) energy density inside the inhomogeneity - a relation we determine explicitly. In combination with continuation arguments and Sturm-Liouville theory, we determine the stability of all constructed pinned fluxons. It follows that if a given system is able to sustain at least one pinned fluxon, a microresistor has exactly one pinned fluxon, i.e. the system selects one unique pinned stable pinned configuration, and a microresonator has at least one stable pinned configuration. Moreover, it is shown that both for micro-resistors and micro-resonators this stable pinned configuration may be non-monotonic - something which is not possible in the homogeneous case. Finally, it is shown that results in the literature on localised inhomogeneities can be recovered as limits of our results on micro-resonators. © 2011 Cambridge University Press.

The stability properties of line solitary wave solutions of the (2+1)-dimensional Boussinesq equation with respect to transverse perturbations and their consequences are considered. A geometric condition arising from a multisymplectic formulation of this equation gives an explicit relation between the parameters for transverse instability when the transverse wave number is small. The Evans function is then computed explicitly, giving the eigenvalues for the transverse instability for all transverse wave numbers. To determine the nonlinear and long-time implications of the transverse instability, numerical simulations are performed using pseudospectral discretization. The numerics confirm the analytic results, and in all cases studied, the transverse instability leads to collapse.

We consider a Josephson junction system installed with a finite length inhomogeneity, either of microresistor or of microresonator type. The system can be modelled by a sine-Gordon equation with a piecewise-constant function to represent the varying Josephson tunneling critical current. The existence of pinned fluxons depends on the length of the inhomogeneity, the variation in the Josephson tunneling critical current and the applied bias current. We establish that a system may either not be able to sustain a pinned fluxon, or - for instance by varying the length of the inhomogeneity - may exhibit various different types of pinned fluxons. Our stability analysis shows that changes of stability can only occur at critical points of the length of the inhomogeneity as a function of the (Hamiltonian) energy density inside the inhomogeneity - a relation we determine explicitly. In combination with continuation arguments and Sturm-Liouville theory, we determine the stability of all constructed pinned fluxons. It follows that if a given system is able to sustain at least one pinned fluxon, there is exactly one stable pinned fluxon, i.e. the system selects one unique stable pinned configuration. Moreover, it is shown that both for microresistors and microresonators this stable pinned configuration may be non-monotonic - something which is not possible in the homogeneous case. Finally, it is shown that results in the literature on localised inhomogeneities can be recovered as limits of our results on microresonators.

The spectral problem associated with the linearization about solitary waves of spinor systems or optical coupled mode equations supporting gap solitons is formulated in terms of the Evans function, a complex analytic function whose zeros correspond to eigenvalues. These problems may exhibit oscillatory instabilities where eigenvalues detach from the edges of the continuous spectrum, so called edge bifurcations. A numerical framework, based on a fast robust shooting algorithm using exterior algebra is described. The complete algorithm is robust in the sense that it does not produce spurious unstable eigenvalues. The algorithm allows to locate exactly where the unstable discrete eigenvalues detach from the continuous spectrum. Moreover, the algorithm allows for stable shooting along multi-dimensional stable and unstable manifolds. The method is illustrated by computing the stability and instability of gap solitary waves of a coupled mode model.