Termination points and homoclinic glueing for a class of inhomogeneous nonlinear ordinary differential equations

Solutions $u(x)$ to the class of inhomogeneous nonlinear ordinary differential equations taking the form \[u'' + u^2 = \alpha f(x) \] for parameter $\alpha$ are studied. The problem is defined on the $x$ line with decay of both the solution $u(x)$ and the imposed forcing $f(x)$ as $|x| \to \infty $. The rate of decay of $f(x)$ is important and has a strong influence on the structure of the solution space. Three particular forcings are examined primarily: a rectilinear top-hat, a Gaussian, and a Lorentzian, the latter two exhibiting exponential and algebraic decay, respectively, for large $x$. The problem for the top hat can be solved exactly, but for the Gaussian and the Lorentzian it must be computed numerically in general. Calculations suggest that an infinite number of solution branches exist in each case. For the top-hat and the Gaussian the solution branches terminate at a discrete set of $\alpha$ values starting from zero. A general asymptotic description of the solutions near to a termination point is constructed that also provides information on the existence of local fold behaviour. The solution branches for the Lorentzian forcing do not terminate in general. For large $\alpha$ the asymptotic analysis of Keeler, Binder $\&$ Blyth (2018"On the critical free-surface flow over localised topography", J. Fluid Mech., 832, 73-96) is extended to describe the behaviour on any given solution branch using a method for glueing homoclinic connections.


Introduction
We investigate solutions to the nonlinear ordinary differential equation, It is assumed that f(0) = 1 and that f → 0 as x → ∞. The rate of decay for large x is a delicate issue and has subtle and important implications for the solution. To highlight this feature of the problem, three particular forcing functions will be examined primarily: a top hat with compact support, a Gaussian and a Lorentzian, given by respectively, where H(x) is the Heaviside function and L is the half-width of the top hat. The Lorentzian function is also known as the witch of Agnesi or as a Cauchy distribution. Assuming that f(x) = f(−x), as is the case for all of the forcings in (1.3), the boundary condition (1.2) may be viewed as providing an even solution over the entire x line, and occasionally it will be helpful to discuss the problem in this context to illuminate some of the key features. Integrating (1.1) directly it is easily seen that provides a necessary condition for a non-trivial solution to exist. For all three forcings in (1.3) the integrand in (1.4) is non-negative and hence non-trivial solutions can only exist for α > 0. The problem is motivated by the study of free-surface flow of an inviscid, irrotational fluid over bottom topography. The forcing function f(x) represents the negative of the topography so that the forcings in (1.3) all correspond to a localised depression on an otherwise flat bottom. In the weakly-nonlinear limit of small forcing, the disturbance to the free-surface induced by the localised topography is described by the forced Korteweg-de Vries equation. The displacement of the free surface from its mean level is given by u(x) governed by (1.1) assuming that the flow is steady and that the speed of the fluid far upstream of the depression is equal to the speed of small amplitude linear waves over a flat bottom (so that the Froude number for the flow is equal to unity).
Recently, Keeler et al [8] considered this problem for the Gaussian forcing. They made a number of observations about the solution space for u(x) that require further mathematical explanation. In particular they presented numerical evidence that there exists an infinite number of distinct solution branches. To place the current work in context, figure 1 shows part of the solution space uncovered by Keeler et al [8], 4 using u(0) to characterise the solutions over a range of values of α. In [8] a traditional boundary layer analysis was used to construct asymptotic approximations both for small α and for large α that approximate the solutions on branch B 0 (in both limits) and on branch B 1 for large α. The rest of the branches, labelled B n for integer n, are not captured by [8]'s asymptotics and this provides one motivation for the present study. The present taxonomy for the solution branches differs from that used in [8] and is motivated by the observation that the solution profiles on branch B n have n local maxima. Since the solution spaces for all three of the forcings in (1.3) share similar qualitative features (but with some key differences), the same taxonomy for the solution branches will be used in each case. Keeler et al [8] provided solid but not conclusive numerical evidence that branch B 1 terminates at its leftmost end at a finite value of α. The present work provides a deeper analysis of this issue.
The layout of the paper is as follows. In section 2 the top-hat forcing is considered; this problem has many of the important features also found for the smooth forcings but with the advantage that the solution can be found exactly. Next in section 3 the importance of the farfield decay rate for a smooth forcing is discussed, and the method for obtaining numerical solutions is described in section 4. In section 5 an asymptotic analysis is presented that supports the termination of the branches B 1 , B 2 etc at finite α and an analysis that indicates that the branch B 0 terminates at α = 0. In section 6 the case of a Lorentzian forcing is examined. Finally in section 7 the method of homoclinic glueing is used to show how the large α solutions can be constructed for a general forcing with a local maximum. The appendices contain further details of the calculations for the homoclinic glueing, a Stokes line analysis for the Lorentzian forcing, and a discussion of a marginal case f(x) = 1/(1 + x 4 ).

Top hat forcing
The top hat forcing, which takes the form given in (1.3), provides an instructive model for the more technically challenging cases (the Gaussian and the Lorentzian forcings), not least because the solution can be obtained exactly in closed form.
A straightforward phase plane analysis nicely illustrates how the key features of the solution space emerge (see Binder [1] for a review of this technique applied to the KdV equation). The unforced phase plane, labelled Σ 1 , corresponds to the homogeneous form of (1.1) and is relevant outside of the top-hat's support where |x| > L. It has a degenerate node at the origin, indicated in figures 2(a)-(c) by an empty circle, with a stable manifold and an unstable manifold on which holds and that are shown each with a broken line. The forced phase plane, labelled Σ 2 , is relevant inside the top-hat support where |x| < L. It has a saddle point at (u, du/dx) = (−α 1/2 , 0) and a centre at (α 1/2 , 0), both of which are shown in figures 2(a)-(c) with filled circles. (Note that the phase portraits Σ 1 and Σ 2 are presented on the same scale.) Trajectories in Σ 2 satisfy for constant c. These are shown with thin solid lines for different c and comprise periodic orbits around the centre enclosed by a homoclinic orbit that connects the saddle to itself. Solutions that satisfy the boundary conditions (1.2) are indicated by thick solid lines in figures 2(a)-(c). In each case, starting from x = −∞ the solution exits the origin and follows the unstable manifold in Σ 1 until x = −L where it jumps instantaneously onto a periodic orbit in Σ 2 . The trajectory jumps instantaneously back onto Σ 1 at x = L and subsequently follows the stable manifold back into the origin as x → ∞. Thus the solutions are smooth everywhere except at x = ±L where the second derivative of u is discontinuous. Various possibilities arise while the trajectory is in Σ 2 , depending on the value of α. A solution that is negative-definite in u can be constructed for any α > 0 by making only a partial excursion along the periodic orbit in the left-half plane of Σ 2 , as is illustrated in figure 2(a). Alternatively a trajectory may execute one cycle of the periodic orbit followed in general by a brief overshoot to make the connection back onto Σ 1 , as is shown in figure 2(b); however this is only possible if the top-hat is sufficiently wide and hence such a solution exists only when α > α * 1 , where α * 1 can be determined precisely and is given below. A countably infinite number of further options arises when α exceeds an increasing sequence of critical values, α = α * n , that can also be written down exactly. For each n the solution executes n cycles of a periodic orbit in Σ 2 followed by an overshoot to connect back to Σ 1 . At the critical values themselves the solution executes exactly n cycles along a periodic orbit in Σ 2 , entering and leaving this plane from Σ 1 at the origin; in Σ 1 itself the solution is given by u = 0 for all x > |L|. This critical case is illustrated in figure 2(c).
The first three solution branches are shown in figure 2(d) together with some sample solution profiles. In all cases the phase plane trajectories are bounded within the homoclinic orbit in Σ 2 ; it follows that −α 1/2 u(x) < 0 on branch B 0 and −α 1/2 u(x) 2α 1/2 on branches B n for n = 1, 2, · · ·. On the periodic orbit in Σ 2 for the critical case, where cn is a Jacobi Elliptic function. Comparing the period of this form to the width of the top-hat we find that for n = 1, 2, · · ·, where K is the complete elliptic integral of the first kind.

Far-field decay for smooth forcings
The unforced, homogeneous form of (1.1) has the general solution that decays at infinity, for arbitrary constant x 0 . Assuming that the generic far-field behaviour of the solution is having a single degree of freedom, namely x 0 in (3.1), which is effectively determined via the choice of α. The large x balance between the first term on the left-hand side of (1.1) and the forcing on the right-hand side, is then also possible but will occur only for certain special values of the parameter, α * n , the behaviour (3.5) involving zero degrees of freedom (in integrating (3.4) to obtain (3.5) both constants of integration must be fixed to ensure the far-field behaviour (1.2) is satisfied). Such a balance neglects the nonlinear term in (1.1) and this is justified provided that (3.2) holds. If (3.2) fails the generic far-field balance is between the nonlinear term and the forcing, given by where we have adopted the negative square root. The positive square root can be excluded on noting that the linearisation at infinity u = ±α 1/2 f 1/2 (x) + U(x) yields If the positive square root is selected then a standard WKBJ analysis of (3.7) yields the two linearly independent unbounded solutions To exclude both of these requires two boundary conditions to be applied at infinity, but this leaves no freedom to enforce the symmetry condition at x = 0 in (1.2). On the contrary, if the negative square root is selected, a single degree of freedom is retained since it is only necessary to exclude the exponentially growing solution to (3.7). The balance (3.6) occurs for the Lorentzian forcing, the third option in the list (1.3), and this case will be examined in section 6.

Numerical computation
Numerical computations for the Gaussian forcing were carried out in [8]. For any of the forcings in (1.3) it is expedient to first rewrite (1.1) as a first order system and then to solve the initial value problem where u = (u, du/dx) T and F = (du/dx, α f − u 2 ) T for some u 0 to be found such that u → 0 as x → ∞ to fulfill (1.2). Thus a solution trajectory in the (u, du/dx) phase plane must ultimately enter the origin and, disregarding the degenerate behaviour (3.4), it must do so in the second quadrant. For a Gaussian forcing, according to (3.3) it will enter the origin along the stable manifold of the degenerate node in the unforced phase plane Σ 1 defined in section 2. The computations for the Lorentzian forcing are particularly challenging as linearising about the far-field decay (3.6) by writing one solution of which, where I 1 is a modified Bessel function, grows exponentially for large x. Hence on shooting from x = 0, any deviation from the required solution will rapidly grow. In computational practice on a finite precision machine any choice for u 0 will result in 'finite-time' blow up with u ∼ u H as x → −x 0 (with x 0 < 0), so that the phase plane trajectory converges to the unstable manifold in Σ 1 . Nevertheless a solution can be detected to good accuracy by using a bisection approach. This is illustrated in figures 3 and 4 for the Gaussian and  Lorentzian forcings respectively. The trajectories were computed by integrating (4.1) using the fourth-order Runge-Kutta method starting in each figure from two carefully selected positive values of u 0 (only a close-up near to the origin is shown). Since in both figures the two trajectories veer either side of the origin, assuming that U(x) depends continuously on u 0 there must exist a u 0 such that the corresponding trajectory reaches the origin and the far-field condition is satisfied.
With confidence that a solution exists, in numerical practice it is more convenient to work with the original boundary value problem statement (1.1) and (1.2). This approach expedites the computation of bifurcation diagrams, for example. To overcome the difficulties presented by the infinite domain and by the algebraic decay of the solution in the far-field we make a change of independent variable. Writing x = cot ζ, with ζ ∈ (0, π), we seek a solution in the form of a Fourier series (e.g. [3, section 17.9]), The coefficients b n are found numerically upon truncation of the infinite sum using a collocation method and employing Newton iterations to solve the nonlinear algebraic equations that arise at the collocation points. The form (4.5) is such that the first boundary condition in (1.2) is automatically satisfied and that u = O(1/x) as |x| → ∞. Therefore with this approach the difficulty associated with (4.4) for the Lorentzian does not arise and accurate solutions that follow the stable manifold in to the origin of the (u, du/dx) phase plane can be computed. This approach also works for Gaussian forcing for which the far-field decay of the solution is like 1/x 2 . Numerical calculations reveal that on a solution branch with a termination point the generic behaviour (3.1) is found at all points along the branch except at the termination point where the singular far-field decay (3.4) is found. This is what was found, for example, in [8] for a Gaussian forcing where at the termination point the decay is superexponential, corresponding to (3.4), and given by Such solutions may be viewed as eigenmodes associated with eigenvalues α * n , contrasting from solutions satisfying (3.3) in existing only for discrete values of α and exhibiting the maximal rate of decay as x → ∞. We leave their relevance to applications as an open question.

Branch termination
The numerical calculations of [8] for a Gaussian forcing suggest that branch B 0 terminates at α = 0 and branches B n , n = 0, 1, 2, . . .terminate at some value α n > 0. In this section we present an asymptotic description of the branch termination in both cases.

Termination at α = 0
As was noted in the Introduction non-trivial solutions to the problem (1.1) and (1.2) for the stated class of forcing functions f(x) exist only if α > 0. Therefore the solution branch B 0 that enters the origin in figure 1 cannot pass into the left-half plane. Keeler et al [8] gave an asymptotic description of solutions on this branch for small α. In fact branch B 0 must terminate at the origin. To demonstrate this, it is helpful to recapitulate some of the key details of the small α analysis.
The expansion proceeds as u(x; α) = u(0; α) + αv(x) + · · ·, noting that u(0; α) = O(α 2/3 ) which follows from the matching carried out below. Substituting into (1.1), at leading order we obtain the linearised form suggests the outer scaling on which the nonlinear term is restored, Under this scaling (1.1) becomes Assuming f(x) = o(1/x) as x → ∞ the right-hand side of (5.4) vanishes to leading order (note that this condition on f is also required for the integral in (5.2) to be defined), and the solution subject to V → 0 as X → ∞ is for constant X 0 , as shown in figure 5(b). (Modulus bars have been included in (5.5) to highlight the singular behaviour at X = 0 discussed in detail below.) This behaves as as X → 0. Matching with the solution on the inner scale yields and as α → 0 + , which agrees with the leading order prediction in [8] for a Gaussian forcing. This approximation is shown with the numerical solution in figure 5(a).
The boundary value problem written in the outer scalings (5.3) best illustrates the way in which the branch terminates: the limit case α = 0 in (5.4) then corresponds to a naïve but natural replacement of (5.4) by in the sense that (5.6) and (5.7) imply dV dX (given the nonlinearity of (5.9) this interpretation should of course be treated with considerable caution as will be seen in the next subsection). Thus the limit profile contains a corner and the branch cannot be continued.

Termination at finite α
As was discussed in section 3 the branches B n for n = 1, 2, . . .terminate at the special values α * n . In this subsection a local asymptotic analysis is presented that describes the termination of an individual branch. To this end we write with 0 < 1, where α * n is one of the special values discussed in section 3 at which the balance (3.4) holds and its value must be determined numerically. We therefore make the implicit assumption that the forcing satisfies (3.2). The choice of sign in (5.11) will be discussed below.
Introducing the expansion and substituting into (1.1) we obtain at leading order in , By the definition of α * n , the far-field decay of u 0 satisfies (3.4). At first order we find The latter of these conditions is required for the matching to be described below and implies that The solution that decays in the far-field is where we have included modulus bars to highlight the singular behaviour at X = 0 to be discussed below. Associated with the solution (5.19) is the requirement that κ < 0 if the plus sign in (5.11) is used, meaning that near to the termination point the branch is such that α > α * n , and the requirement that κ > 0 if the minus sign in (5.11) is used so that α < α * n local to the termination point. These requirements ensure a match with the solution on the inner scale. It follows that a sufficient condition for the existence of a fold in the solution branch is that κ > 0. We may infer from the numerical results shown in figure 1 that κ > 0 for the Gaussian forcing. It should be emphasised that while the present analysis gives a self-consistent description of the behaviour close to a termination point, it does not prove their existence even for forcings that satisfy the far-field decay condition (3.2). It may preclude the existence of a termination point, however, if (3.2) is not satisfied. Key to the latter remark is the existence of two possible large x balances for forcings that satisfy (3.2), namely (3.3) and (3.5), on which the analysis presented in this subsection depends. Forcings that do not satisfy (3.2) have only one possible large x balance as discussed in section 3.
Since u(0) ∼ u 0 (0) + u 1 (0) + · · · the branches in figure 1 are locally linear and enter the termination points with finite slope. In common with the small α case discussed in section 5.1 the limit profile at each termination point, on the outer scale (5.19), has a corner at X = 0 with the jump in slope dv 0 dX Notably, and in contrast to (5.10), the jump is not given simply in terms of the mass of f(x) but relies on the numerically determined constant κ. This underscores the danger, alluded to in the previous section, of a naïve replacement of the right-hand side of the problem on the outer scale, here (5.18), with Mδ(X), where M is the mass of the forcing given in (5.2), since the jump in slope at X = 0 depends on the solution of the problem on the inner scale according to (5.20). The preceding remarks can be placed on a firmer footing by noting that on taking the limit α → 0 in (5.4) the right-hand side formally approaches a delta function (e.g. Stakgold and Holst [9], theorem 2.2.4). On the contrary, writing (1.1) for the outer scaling results in (5.18) with on the right-hand side, and this does not approach a delta function in the limit → 0. Consequently, branch B 0 can be described by naïvely replacing the right-hand side of (1.1) with Mδ(x) and then following the type of phase plane analysis reviewed by Binder [1], but the remaining branches B n for n 1 cannot be described in this way.

Lorentzian forcing
The numerically computed solution space for the Lorentzian forcing in (1.3) is similar in structure to that found for the Gaussian (see figure 1) with the crucial difference that the higher order branches B 1 , B 2 , etc do not terminate at finite α. As for the Gaussian forcing there is a B 0 branch of negative-definite solutions that exists for all positive α and which terminates at α = 0 as described in section 5.1. We do not expect to find branches that terminate at non-zero α for reasons discussed in section 3. In fact we find that the branches B 1 and B 2 are in this case connected continuously as can be seen in figure 6. The insets in this figure show typical solution profiles some way along the upper and lower parts of the branch. Note that although the comment at the end of section 3 leaves open the possibility that eigenmode solutions exist for special values of α (as in the Gaussian case) corresponding to the choice of the plus sign in (3.7), our numerical computations suggest that such solutions do not, in fact, exist.
Following the success of the boundary-layer analysis of [8] on branch B 0 for the Gaussian forcing, we are motivated to attempt a similar large α analysis for the Lorentzian, and this is considered in the following subsection. As for the Gaussian such an analysis cannot capture the higher order branches. These will be discussed in section 7.

Asymptotic approximation when α 1
For large α we rescale by writing u = α 1/2 W(x) so that (1.1) becomes where μ = α −1/2 . We seek an asymptotic expansion in the form Substituting into (6.1) and working at successive orders, we obtain In general for n 1, and it is straightforward to show that W n (x) = p 2n (x)/(1 + x 2 ) (1+3n)/2 , where p 2n (x) is a polynomial of degree 2n. For any n an important observation is that W n satisfies the boundary conditions (1.2) irrespective of the choice of sign in (6.3). This is reminiscent of the problems discussed by Chapman et al [5] and references therein, whereby for some ordinary differential equations with a small parameter, a simple asymptotic solution can be constructed that satisfies the required boundary conditions at any order and for all x when no such solution to the problem in fact exists. The difficulty can be traced to the fact that the expansion disorders in the neighbourhood of the singularities in the leading order term in the asymptotic expansion (W 0 here) extended into the complex plane. If a Stokes line emanating from one of these singularities crosses the real axis, an exponentially small 'beyond all orders' term is in general switched on at the point of crossing and this term eventually grows to corrupt the original expansion. In appendix B we provide the relevant Stokes line analysis for the present problem. The conclusion is that if the plus sign in (6.3) is chosen then the expansion (6.2) is corrupted in the manner described. For the minus sign no such difficulty arises and the asymptotic expansion (6.2) is valid for all x.
Selecting the minus sign for W 0 in (6.2) we find u(0) ∼ −α 1/2 + 1/2. This approximation is shown with a broken line in figure 6. As this figure suggests, further possibilities arise for the large α asymptotics in which (6.2) is regarded as an outer expansion to be matched to an inner boundary-layer solution around x = 0 that describes a cluster of localised waves. These localised waves may be viewed as the connection of individual solitary-wave structures that are each described by a homoclinic orbit in an appropriately defined phase space, as will be discussed in the following section.

Homoclinic glueing
In this section we aim to describe asymptotic forms that approximate the solutions in the limit of large α for a fairly general class of forcing functions f. It will be convenient to think in terms of symmetric solutions to (1.1) that are defined on the whole of the real line and that decay as |x| → ±∞. To motivate the construction we rewrite (1.1) on a boundary-layer scale by introducing the new variable x = α −1/4 y and setting u(x) = α 1/2 U(y) to obtain where δ = α −1/4 1. Expanding the solution as U = U 0 (y) + δ 2 U 1 (y) + · · ·, and replacing the right-hand side by its Taylor expansion f(δy) = f(0) + (1/2)δ 2 f (0)y 2 + O(δ 3 ), at leading order we find (since f(0) = 1) This has three bounded solutions of interest, that correspond to equilibria in the (U 0 , dU 0 /dy) phase plane (i) and (ii) and a homoclinic orbit connecting (−1, 0) to itself in the same plane (iii). It is symmetric about y = 0 and has the property as y → ∞. (7.4) Its graph in physical space has a classical solitary-wave shape (see, for example, Billingham and King [2]), and this suggests representing the wave-like parts of the B n branch solutions by a collection of these homoclinics. The strategy is as follows: for n odd, we seek to glue together n of the homoclinics (iii) in (7.3) via asymptotic matching; for n even, n − 1 homoclinics are glued together either side of a central region in which solution (ii) in (7.3) predominates. We consider odd and even numbers of homoclinics separately in the following subsections. In both cases the analysis is predicated on the assumption that f (0) < 0, so that the forcing has a local maximum at the origin (or, by a suitable shift, at any x location). This condition is fulfilled both by the Gaussian and the Lorentzian forcings.

Odd number of homoclinics
Continuing with the analysis, at next order we have Henceforth in this subsection it is assumed that U 0 is given by the homoclinic (iii) in (7.3). The general solution to (7.5) is given in appendix A in terms of a symmetric and an antisymmetric complementary function Φ s and Φ a , and particular integral form Φ (2) p that satisfies Φ (2) p (0) = dΦ (2) p /dy(0) = 0. The solution that satisfies the boundary conditions U 1 (0) = U 10 , dU 1 dy (0) = 0 (7.6) is given by where U 10 is a constant that will be determined later. We note from (A.8) that For a single homoclinic we require to exclude the exponential growth as y → +∞. This is the case considered in [8] for the Gaussian forcing. The solution for U(y) is then matched to the solution on the outer scale, where x = O(1). 5 More generally, using the asymptotic forms given in appendix A we have as the required matching condition 2y + · · · + · · · (7.10) as y → +∞, where According to (7.9) the single homoclinic has Λ 1 = 0. If Λ 1 > 0 then a second homoclinic is initiated at 12) whereupon (7.10) becomes 2y 1 + · · · + · · · (7.13) as y 1 → −∞. A third homoclinic is initiated in y < 0 to maintain symmetry.

Remark 1.
The shift Y 1 (δ) has been chosen to make the correction term δ 2 Λ 1 e √ 2y in (7.10) of O(1) and so that the two exponential terms in (7.10) are effectively interchanged in (7.13) to effect the matching between the homoclinics.
For y 1 = O(1) we have the expansion U = U H 0 (y 1 ) + δ 2 θ 1 (y 1 ; log(1/δ)) + · · · , (7.14) where U H 0 was given in (7.3) and provided that δ 2 Y 1 1, a restriction that will be discussed in more detail below. Henceforth the subscripts on θ will label the homoclinic sequence rather than the δ expansion. Also, as is suggested by the notation, here and subsequently we shall lump all additional logarithmic factors into θ i in order to obtain algebraic rather than simply logarithmic accuracy.

Remark 2.
The inner expansion at O(δ 4 ) will lead to a δ 2 e √ 2y 1 term in (7.13). This will simply trigger a dU H 0 /dy 1 complementary function in the solution to (7.15), which corresponds to an O(δ 2 ) translation in y 1 . This can be safely ignored since we do not seek to determine O(δ 2 ) corrections to the locations of the maxima. Inspecting (7.13) and (7.15) we decompose θ 1 as where φ and the ψ k , k = 0, 1, 2 satisfy the problems stated in appendix A. To perform the matching we demand that The given stipulations for the · · · in both (7.17) and (7.18) are made to remove the translational invariance alluded to above to ensure a unique solution.

Even number of homoclinics
The first pair of homoclinics is located at y = ±Y 1 (δ), where Y 1 is to be found (note that Y 1 now differs from that given in section 7.1). Sufficiently close to y = 0 the expansion U = −1 + δ 2 U 1 (y) + · · · (7.27) holds, where the leading order term corresponds to (ii) in (7.3). Note that the choice (i) in (7.3) was ruled out in [8] and is ruled out here for the same reason. At first order with solution where the constants of integration have been set to ensure a match with the homoclinics at y = ±Y 1 .
Remark 3. (7.29) appears to be inconsistent with the expansion (7.27); in fact Y 1 will be such that the 1/δ 2 terms in (7.29) are of O(1/δ) approximately, a statement that will be made more precise below, so |δ 2 U 1 | = O(δ) 1. Table 1. Comparison between numerical and asymptotic results for the homoclinic glueing. In the table u 0 is the numerical estimate and u A 0 is given by the pertinent asymptotic formula (7.39). The location of the first homoclinic maximum is given numerically by Y H and the asymptotic estimate Y 1 is taken to be the numerical solution of (7.12) (with Λ 1 given by the equation in (7.20)) for the triple homoclinic and as the numerical solution of (7.33) for the double homoclinic. The relative errors are δu = |(u 0 − u A 0 )/u 0 | and Setting y = Y 1 + y 1 , and inserting (7.29), (7.27) becomes 2Y 1 e − √ 2y 1 + · · · , (7.30) which motivates writing u = U H 0 (y 1 ) + δ 2 θ 1 (y 1 ; log(1/δ)) for y 1 = O(1), where ψ 0 (y 1 ) + 2Y 1 ψ 1 (y 1 ) + ψ 2 (y 1 ) (7.31) (the functions φ and ψ k , for k = 0, 1, 2, were defined in section 7.1). Once Y 1 is determined below, it can be confirmed a posteriori that each of the terms in (7.31) are at worst logarithmic in δ. Using the details given in appendix A, as y 1 → +∞.
The double homoclinic corresponds to Λ 1 = 0 or 12 2), shown with a solid line. The circles correspond to u = α 1/2 U 0 + U 1 with U 0 given by U H 0 or −1, and U 1 given by (7.7) or (7.29) for the triple/double homoclinic respectively. The crosses correspond to u = α 1/2 U 0 + θ 1 with U 0 given by U H 0 and θ 1 given by (7.16) or (7.31) for the triple/double homoclinic respectively. For the triple and double homoclinics Y 1 was taken to be the numerical solution of (7.12) [with Λ 1 given by the equation in (7.20)] and (7.33), respectively. in which case Otherwise the analysis continues in a manner similar to that presented in section 7.1. In this case the sequence is established as andŷ n = Y n+1 − Y n +ŷ n+1 for n = 1, 2, · · ·. The sequence for each successive n corresponds to the value Y 1 such that Λ n = 0 and has homoclinics at y = ±Y k for k = 1, . . . , n. We find Taking into consideration remark 3 we may now check the validity of the expansion (7.27). By taking the square root of (7.33) it is clear that the 1/δ 2 terms in (7.29) are in fact of O(λ) where λ = 1 δ log 1/2 1/δ . (7.37) The expansion (7.27) should therefore be adjusted accordingly; this adjustment is permitted due to the linearity of the perturbation problems at each order of approximation. The odd/even homoclinic analysis above is valid provided that the number of homoclinics n is such that δ 2 Y n 1. Given (7.26) and (7.36) this implies n 1 δ 2 log(1/δ) (7.38) as the condition that all the homoclinics are located where |x| 1. Since f (0) = −2 for both the Gaussian and the Lorentzian, the results in the previous subsections give the following approximations which may be applied to either case, In table 1 we compare these asymptotic predictions with numerical calculations for the triple homoclinic for the Gaussian and the Lorentzian, and the double homoclinic for the Gaussian. In figure 7 we show a comparison between the asymptotic homoclinic glueing predictions and numerical solutions for the Gaussian forcing that demonstrate strong agreement between the two.

Discussion
We have analysed solutions to the problem (1.1) and (1.2) for the case of a top hat forcing, a Gaussian forcing and a Lorentzian forcing, with particular attention paid to the limits of small and large α. We have presented an asymptotic construction to provide supporting evidence for the existence of termination points on the solution branches for forcing functions which decay in the far-field faster than 1/x 4 , which includes the Gaussian forcing. We have also presented an asymptotic description of the large α solution profiles using the method of homoclinic glueing, which can be applied to any smooth forcing with a local maximum. The structure of the solution space is similar for all three forcing functions, each with an apparently infinite number of solution branches with qualitatively similar features in the solution profiles on each branch. Of particular note, however, is the presence of termination points for the Gaussian forcing on all branches, and the linking together of branches B n , B n+1 for the Lorentzian forcing (which decays more slowly that 1/x 4 in the far-field). Some further insight into these different characteristic features can be obtained by attempting to continuously deform one forcing function into another. This can be achieved by considering the hybrid forcing function for 0 a 1. As was noted above, the generic far-field behaviour for u(x) in (1.1) corresponds to blow-up at the finite value x = −x 0 via (3.1). Figure 8 shows contours of the blow-up point x 0 < 0 for the forcing (8.1) at three sample values of a. Solution branches can be discerned on which formally x 0 = −∞ to satisfy the far-field condition (1.2). Of particular note when a < 1 are the two saddle points in the contour map that are located for a = 0.7 roughly at (α, u 0 ) = (11, 3.6) and (48, 7.5). These saddle points persist on decreasing a to zero, and indeed the contour plot for the Gaussian forcing is qualitatively similar to figure 8. As a is increased towards unity the tips of the two solution branches move toward each other. They pinch together at the rightmost of the two saddle points when a = 1 to form the continuous solution branch labelled B 1 , B 2 for the Lorentzian forcing in figure 6. As has been discussed, the existence of termination points at non-zero α appears to hinge on the decay rate of the forcing in the far-field relative to the inverse fourth power (see equation (3.2)). A forcing function for which f(x) = O(1/x 4 ) as x → ∞, for example, presents a marginal case. For this forcing, the large x behaviour of the solution is constituting a balance between all three terms in (1.1). The branch B 1 for the forcing (8.2) is shown in figure 9. The inset suggests that the branch terminates at α ≈ 19.9; in fact the branch spirals inwards toward the termination point beyond the last point reached by our numerics (see appendix C for details). This underscores the subtle behaviour that can be found in problems of this type. While we have used the particular class of equation (1.1) as an example to illustrate the idea of a termination point and to demonstrate the use of the method of homoclinic glueing, we believe that this class is simple enough to act as a paradigm for a much broader set of problems. Finally we note that the homoclinic glueing analysis presented here is valid provided that the number of homoclinics satisfies condition (7.38). If this condition is violated then a different approach is needed. This is the subject of ongoing work.
can be readily obtained readily by making the substitution t = tanh(y/ √ 2). The solution is found to be For reference we note the asymptotic properties of these functions. Herein we provide details of the first order problems satisfied by the functions φ(y 1 ) and ψ k (y 1 ) for k = 0, 1, 2 which appear in the general solution for θ 1 in (7.16). The problem for φ is Using the results from above, the general solution may be written as for constants C a , C s . Inspecting the asymptotic forms (A.8) and (A.9) we see that to fulfill the glueing condition (7.17) we must set C a = 0 and C s = −16 and so the solution for φ is even in y 1 .
The functions ψ k satisfy the problems for k = 0, 1, 2. The solutions are: where the D (k) a and D (k) s are arbitrary constants and We note the following asymptotic properties. As y 1 → +∞ and as y 1 → −∞ ,

Appendix B. Stokes line analysis for the Lorentzian forcing
The difficulty noted by Chapman et al [5], for example, requires an analysis of the Stokes lines that emanate from the singularities of the leading order term in the expansion (6.2) extended into the complex plane. With this in mind, when referring to the analysis in section 6.1 we shall replace x with z ∈ C. The expansion (6.2) is a divergent asymptotic series whose terms are generated by the recurrence relation (6.4). Following Chapman et al [5] we optimally truncate the asymptotic series at its smallest term, writing where R N (z) is a remainder term and m = ±1 corresponding to the choice of sign made in (6.3). The optimal truncation level N follows from knowledge of the large n behaviour of W n (z). We make the usual ansatz, writing (see Dingle [6]), as n → ∞, where Γ(z) is the Gamma function and the functions A(z), γ(z) and the singulant σ(z) are all to be found. Substituting (B.2) into (6.4) the balances at leading order, first order and second order determine that (see Keeler [7]) where means d/dz, and that γ(z) = −1/6 and A(z) = Λ/(σ ) 1/2 , where Λ is the Stokes multiplier that will be determined below. Since W 0 has singularities at z = ±i, then W n (z) will also have singularities at these locations, for all n. We shall focus on the singularity at z = i, the analysis for z = −i being similar. Integrating (B.3), the singulant takes the form where the lower integration limit has been chosen so that σ(i) = 0. Finally, consistency of W n (z) for large n between the forms given in (B.2) and in (B.1) demands that (Keeler [7]) Λ = 2 7/6 π 1/2 3 2/3 Γ Stokes lines emerge from the points z = ±i where σ vanishes and hence, according to (B.2), the late form of W n is singular. According to Dingle [6] the Stokes lines are traced by delineating the curves in the complex plane on which successive terms of the late order form (B.2), namely W n and W n+1 , have the same phase. Equivalently (see, e.g. [5]) on a Stokes line, Im(σ) = 0. (B.6) The angle at which the lines emerge from the singularities is determined as follows. We note that σ ∼ 2 9/4 3 i 7/4 (z − i) 3/4 (for m = −1), σ ∼ 2 9/4 3 i 11/4 (z − i) 3/4 (for m = 1), (B.7) as z → i. Since the original problem posed on the real line is symmetric about x = 0, it is natural to take the branch cuts that stem from the branch points at z = ±i to extend up the imaginary axis from z = i and down the imaginary axis from z = −i respectively. If local to z = i we write z − i = Re iψ and σ = re iθ , then we should insist that The numerically computed Stokes lines in the upper half plane are shown in figure B1(a) together with the asymptotic approximation (B.9). We note that the thickness of the Stokes layers about each Stokes line is of order μ 1/4 |σ| 1/2 . Since |σ| grows like |z| 1/2 for large x, the Stokes layer thickness grows as x 1/4 , and hence the Stokes layer cannot impinge on the real line. It follows that for m = −1 the Stokes phenomenon can be ignored and the expansion (6.2) holds for all real x.
The situation is different for m = 1. In this case the local form (B.7) leads to the conclusion that there is a Stokes' line along the imaginary axis between z = i and z = −i. An analysis similar to that in Chapman et al [5] (see Keeler [7] for details) shows that the exponentially small remainder term R N (x) ∼ 2 3/4 πΛ μ −5/12 e −ρ(2/μ) 1/2 (1 + x 2 ) 1/8 cos ψ(x), (B.10) where ψ(x) = 2/μ 1/2 x 0 (1 + p 2 ) −1/4 dp, is activated on the real axis where the Stokes line crosses it at x = 0. This term grows algebraically in x so eventually the expansion (6.2) breaks down and there is therefore no solution of the original boundary value problem which is approximated by (6.2) for all x. with du 1 /dx(0) = 0 and u 1 (x) → 0 as x → ∞. For large x, we have u * ∼ ϕ + (ξ)/x 2 and the complementary functions for (C.6) are 1/(x 2 x p ± ). Numerical computations suggest that α * is such that p ± are a complex conjugate pair, so that for large x, u ∼ 1 for complex constants A ± (with A + = A − ) and τ = (8 √ 9 + α − 25) 1/2 . A single relation is needed between the constants A ± to have a boundary value problem for u 1 .
The non-uniformity of the expansion (C.7) implies the presence of an outer region in which x = −2/5 X and u = 4/5 U with X = O(1) and U = O(1). Writing U = Φ/X 2 and Θ = log X, we have that is two equations which provide the relation between A + and A − alluded to above, and a condition to determine Θ 0 , on solving the boundary value problem for u 1 .
According to (C.11) the rescaling → e −5π/τ , arg A ± → arg A ± ± 2π (C.12) leaves the solution for u 1 unchanged. Therefore in this case the solution branch spirals into the termination point at α = α * in the manner sketched in figure C1.