Semi-infinite travelling waves arising in a general reaction–diffusion Stefan model

We examine travelling wave solutions of the reaction–diffusion equation, ∂tu=R(u)+∂xD(u)∂xu , with a Stefan-like condition at the edge of the moving front. With only a few assumptions on R(u) and D(u), a variety of new semi-infinite travelling waves arise in this reaction–diffusion Stefan model. While other reaction–diffusion models can admit semi-infinite travelling waves for a unique wavespeed, we show that semi-infinite travelling waves in the reaction–diffusion Stefan model exist over a range of wavespeeds. Furthermore, we determine the necessary conditions on R(u) and D(u) for which semi-infinite travelling waves exist for all wavespeeds. Using asymptotic analysis in various distinguished limits of the wavespeed, we obtain approximate solutions of these travelling waves, agreeing with numerical simulations with high accuracy.

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the one-dimensional domain x ∈ R at time t > 0. Notably, there has been particular interest in semi-infinite travelling waves, in which u(x, t) is a monotone decreasing function that terminates at u = 0 when at some finite spatial position x(t). These travelling waves are useful in applications where a well-defined 'edge' of a travelling wave is an essential component of the modelling framework [10][11][12]. While many existing reaction-diffusion models have been shown to produce a semi-infinite travelling wave for a single wavespeed, c * [2][3][4][18][19][20][21][22][23], there are only a handful of modelling frameworks giving rise to semi-infinite travelling waves for range of wavespeeds.
One particular reaction-diffusion modelling framework giving rise to semi-infinite travelling waves over a range of wavespeeds is via the incorporation of a moving boundary [8][9][10][11]24], whereby x ∈ (−∞, L(t)] and L(t) evolves based on a Stefan-like condition at the edge of travelling wave. For particular choices of linear [8][9][10]24] and degenerate diffusivities [11], D(u), semi-infinite travelling waves exist for all wavespeeds c ∈ [0, c * ], where the value of the critical wavespeed c * depends on D(u). However, generalisations of the models presented in [10,11] for a broader class of reaction functions, R(u) and nonlinear diffusivities, D (u), has yet to be considered.
In this work, we consider a general reaction-diffusion model with a moving boundary, which we refer to as the reaction-diffusion Stefan (RDS) model. Specifically, we determine the necessary conditions for R(u) and D(u) so that semi-infinite travelling waves exist over a range of wavespeeds. Using asymptotic analysis in the limit where c 1, we obtain approximate solutions of these travelling waves, agreeing with numerical simulations with high accuracy. Additionally, we determine the approximate relationship between c and the Stefan parameter, κ, which relates the concentration flux, −D(u)∂ x u(x, t), with the speed of the moving boundary. Along with the asymptotic analysis performed when c 1, corresponding to when κ 1, we also consider approximate solutions of travelling waves for κ 1. Depending on the kinetics of R(u)D(u) near u = 0, we show that having κ 1 is equivalent to c approaching a finite wavespeed, or that c 1. In particular, we outline the necessary conditions for R(u) and D(u) so that semi-infinite travelling waves exist for all wavespeeds (i.e., c * = ∞). For both c * finite and infinite, we determine asymptotic approximations for both the travelling wave and the κ(c) relationship in the limit where κ 1.

Travelling waves in the reaction-diffusion Stefan model
We consider the following non-dimensional reaction-diffusion model, describing the concentration u ∈ [0, 1], on the spatial domain x ∈ (−∞, L(t)] with a Stefan-like condition at the moving boundary x = L(t): The Stefan condition relates the speed of the moving front, dL/dt, to the concentration flux, −D(u)∂ x u(x, t) via the Stefan parameter, κ. We will also impose that These bounds indicate that degenerate diffusion, i.e. D = 0 or D = ∞, can only occur if u = 0, and that the reaction R(u) can only be a source term. The main solutions of the system (1)-(3) that we will examine are travelling wave solutions. Therefore, we transform the PDE system into travelling wave coordinates via z = x − L 0 − ct, where z ∈ (−∞, 0] and c 0. Noting that when x = L(t), we have that L(t) = L 0 + ct, and hence, dL/dt = c. By denoting the concentration flux as we obtain, by transforming (1) into travelling wave coordinates and multiplying by D(u), Provided that D(u)du/dz → 0 as u → 1 − , the system of first-order nonlinear differential equations (5) and (6) is coupled to a set of boundary conditions arising from (3): We will focus on the heteroclinic trajectory in the (u, Q) phase plane beginning at (1, 0) and terminating at (0, −c/κ) to agree with the boundary conditions listed above. In doing so, we determine not only the relationship between Q and u, but also the relationship between c and κ. Therefore, by dividing (6) by (5), we have that

Some explicit solutions of the heteroclinic trajectory Q(u)
While by no means an exhaustive list, there are some straightforward solutions of (8) for particular choices of R(u)D(u). In doing so, we are also able to obtain the explicit relationship between c and κ. The simplest case to consider is when R(u)D(u) = A, where A > 0. In this case, (8) becomes a separable differential equation whose solution is where W(·) is the Lambert-W function [25]. Therefore, by evaluating this expression at u = 0, we retrieve the relationship between c and κ: However, since A > 0, we have, from (8), that dQ(u)/du → −∞ as u → 1 − . Furthermore, since D(1) is positive and finite, the travelling wave u(z) will evolve on some finite interval [Z 0 , 0] and u(z) ≡ 1 for z Z 0 . Should a travelling wave be required to evolve on the semi-infinite interval (−∞, 0], we must further impose that R(1) = 0. However, as we will show, having R(1) = 0 is a necessary but not sufficient condition in obtaining travelling waves on the semi-infinite interval (−∞, 0]. A simple reaction-diffusion model that exhibits travelling waves on a semi-infinite domain Therefore, as u → 1 − , du/dz ∝ u − 1, confirming that travelling waves in this reaction-diffusion model decay exponentially to u = 1 as z → −∞.

Approximating the heteroclinic trajectory Q(u) for c 1
For more general R(u)D(u) function forms, an explicit solution of (8) is not always possible. Instead, we consider the solution of (8) first in the limit where 0 < c 1 by expanding Q(u) as a regular perturbation expansion in c, i.e. Q(u) = q 0 (u) + cq 1 (u) + O(c 2 ). Substituting this expansion into (8) provides the following two-term equations: Solving (12) and (13) yields From (4), we have that ξ(u) is non-negative, bounded, and monotone decreasing for u ∈ [0, 1]. Furthermore, evaluating (14) at u = 0 provides a two-term approximation for the wave speed c as a function of the Stefan parameter κ when c 1: Therefore, we can obtain two approximate relationships in the limit where the wavespeed is small: the relationship between the travelling wave and its flux, as well as the relationship between the wavespeed and the Stefan parameter.
By substituting (16) into (14), we obtain where Therefore, in the limit where u → 0 + and c 1, we have, by employing the initial condition Crucially, we note that D(u) must be o(u −1 ) for the right-hand side of (19) to converge, thereby allowing z → 0 as u → 0 + . By approximating D(u) ∼ Au p for u → 0 + , where p > −1, we have that the leading edge of the travelling wave can be determined implicitly as Thus, by knowing the approximate expansions of D(u) and R(u)D(u) near u = 0, we are able to approximate u(z) for u → 0 + .

Approximating the travelling wave u(z) for c 1 and u → 1 −
Having approximated the travelling wave u(z) near its leading edge, i.e. when u 1, we now examine the far-field behaviour of u(z) by considering the limit where u → 1 − . From the conditions imposed in (4), we have that D(1) is strictly positive and finite, as well as We note that in the special case where γ = 1, du/dz ∝ 1 − u at leading order, implying that where Z 0 is an arbitrary constant. Next, when 0 γ < 1, we rearrange (21) and obtain Since γ < 1, the bracketed term in the integrand in (23) is strictly positive for u ∈ [0, 1] and is close to 1 for c 1 and u → 1 − . Therefore, we obtain the implicit far-field approximation Furthermore, we note that when u = 1, z = Z 0 , indicating that when 0 γ < 1, travelling waves will evolve over a finite interval. Finally, we examine the case where γ > 1. From (23), we identify that the integrand will now have two singularities in [0, 1]: Therefore, for the travelling wave to deviate significantly away from u = 1, we must have c ≡ 0 when γ > 1. This key result indicates that in order for travelling waves in the RDS model to have a semi-infinite domain with non-zero wavespeed, then

Comparison between asymptotic approximations and numerical solutions
Now knowing the approximate forms of the travelling wave, u(z), near its leading edge and near u = 1 for small wavespeeds, we compare the numerically-determined travelling waves with these approximations for various reaction-diffusion models. To compute the numerical solution of (5)-(7) for various R(u)D(u) and c inputs, we use ode15s in MATLAB. The numerically-determined travelling waves are computed over a sufficiently large range of z, with initial conditions very close to the far-field boundary condition (7). For numerical solutions shown here, z ∈ [−10 4 , 0], u(−10 4 ) = 1-10 −9 , and Q(−10 4 ) = 10 −9 . The resulting numerical travelling wave is then translated in z to ensure that u(0) = 0. As a first example of common reaction-diffusion combinations examined in the context of travelling waves, we consider the generalised Porous-Fisher (GPF) model [1,3,4,12,[21][22][23], which generalises the common Fisher-Kolmogorov reaction-diffusion model [5,6] to include nonlinear diffusion: From (14), we obtain While an explicit solution for 1 u ξ(s)ds is not possible for general r, employing (20) implies that for u → 0 + , where and can be easily computed numerically. Additionally, as u → 1 − , the GPF model has, from (22), the far-field approximation Using both of these asymptotic approximations, we see in figure 1 good agreement with the numerically-determined travelling wave solution. While the two-term approximation of Q(u) via (14) must be numerically integrated, we also see good agreement between the two-term approximation and the numerically-determined heteroclinic trajectory for small values of c.
To contrast with the GPF model, we now consider a new reaction-diffusion model that exhibits travelling waves over a finite interval: We refer to this reaction-diffusion model as the fast diffusion, slow reaction (FDSR) model, based on the kinetics near u = 0. While we do not attempt to provide a physical justification for the FDSR, the topics of fast degenerate diffusion are often studied in the context of travelling wave solutions [1,3,12]. The exact solution for Q(u) is shown in (9) with A = 1; however, obtaining u(z) for the FDSR model cannot be explicitly obtained. Using (20) and (24), we have that for c 1, indicating that the FDSR travelling waves do indeed evolve over a finite interval. As shown in figure 1, these asymptotic approximations agree well with the numerically-determined travelling wave solution. Furthermore, the approximation of Q(u) via (14) agrees very well with the exact solution shown in (9). For both the GPF and FDSR models, we can also determine the approximate relationship between c and κ when c 1 via (15): As seen in figure 1, these approximations for κ(c) agree well with the numerically computed values of κ.

Approximating reaction-diffusion Stefan travelling waves for κ 1
Now knowing the the approximations of u(z) and Q(u) in the limit where c 1, i.e. κ 1, a natural extension is to examine travelling waves of the RDS model when the Stefan-like parameter κ is large. It is well-reported that beyond a critical wavespeed, c = c * , travelling waves of reaction-diffusion models are no longer semi-infinite, implying that κ → ∞ as c → c * − [10,11,[21][22][23]. In the context of the heteroclinic trajectory Q(u), which solves (8), we have that Q(0) = 0 for c c * . While this critical wavespeed is normally referred to in the literature as the minimum wavespeed for which smooth travelling waves exist [1][2][3][4][18][19][20][21][22][23], we will interpret c = c * as the maximum wavespeed which produce semi-infinite travelling waves. Additionally, the heteroclinic trajectory Q(u) will 'flatten' and decrease in absolute magnitude as c increases, due to the key result that travelling waves become more diffuse and spread-out as the wavespeed increases [1]. We will first show what conditions on R(u)D(u) are required for c * = ∞, as well as the approximations of Q(u) for c 1. Following this analysis, we will examine instances of R(u)D(u) in which c * is finite, as well as approximations of Q(u) and u(z) for c → c * − .
To determine necessary conditions for c * = ∞, we first examine a desingularised version of equations (5) and (6) via the change of variables This change of variables allows us to examine the equilibrium (u, Q) = (0, 0) without the issue of having degenerate D(u) as u → 0 + , while leaving Q(u) unchanged [1,20]. Furthermore, it immediately follows that in order for (u, Q) = (0, 0) to be an equilibrium of the desingularised system, we must have R(0)D(0) = 0. Thus, for R(0)D(0) = 0, no finite value of c produces a heteroclinic trajectory with Q(0) = 0. However, having R(0)D(0) = 0 is a necessary but not sufficient condition for c * to be finite. By examining the Jacobian J of the desingularised system at (u, Q) = (0, 0), we find that However, this local analysis fails when |α| = ∞, implying that Q(0) = 0 and providing a second case for when c * = ∞. To summarise, the critical wavespeed c * = ∞ when at least one of the two following conditions hold:

Approximating the heteroclinic trajectory Q(u) for c 1
We first examine the heteroclinic trajectory Q(u) for c 1, assuming that the conditions shown in (37) hold. By rescaling Q(u) = cφ(u), (8) becomes Since the left-hand side of (38) is o(1), we anticipate that the boundary conditions will not necessarily be satisfied without additional rescaling of variables. Therefore, we first consider the outer solution of (38) by expanding φ(u) hence, the outer solution for Q(u) is There are two problems that can arise in the outer solution. The first is if R = R(1)D(1) > 0, implying that a boundary layer near u = 1 must exist to satisfy the boundary condition φ(1) = 0. By rescaling u = 1 − c −2 v and denoting φ = ψ(v) in the boundary layer, we obtain the leading-order ODE which has the solution ψ The second problem that can arise is due to the conditions that can hold in (37) as u → 0 + . As a result, the outer solution no longer remains a well-ordered asymptotic expansion, which is resolved by rescaling (38) as u → 0 + . By assuming R(u)D(u) ∼ Bu q for u → 0 + , where q ∈ (−1, 1), a balance of all terms in (38) is achieved via the rescalings which yields the leading-order ODE in the boundary layer near u = 0: A special case of (43) is when q = 0, where the solution is Φ(χ) ≡ −B and Q(u) therefore does not change in this boundary layer. While (43) does not have an explicit solution for general q ∈ (−1, 1), we are still able to provide some insight about the composite leading-order approximation of Q(u) for large c 1: In particular, we have that which provides a leader-order power law relationship between κ and c for large wavespeeds.

Approximating the heteroclinic trajectory Q(u) for c → c * −
As shown in (37) It immediately follows that if α > 0, (u, Q) = (0, 0) changes from a stable spiral point to a saddle point at c * = 2 √ α, thereby providing an explicit expression for the critical wavespeed [1,20,21,24]. As we impose that R(u)D(u) 0, we do not consider the case where α < 0.
For α = 0, the eigenvalues of J (0,0) are (0, −c * ) for all positive wavespeeds, resulting in a fixed point with degenerate stability. Therefore, we cannot use linear stability analysis to determine c * and instead define the critical wavespeed as the minimum value of c that solves (8) and has Q(0) = 0. Additionally, from the eigenvalues of J (0,0) , we have that Q (0) = −c * at the critical wavespeed. While neither Q(u) or c * can be explicitly determined for all choices of R(u)D(u) with c * finite, we are nevertheless interested in examining both Q(u) and κ(c) as c → c * − . To do this, we set c = c * − ε, where ε 1, whereby (8) becomes By expanding Q(u) as a regular perturbation expansion in ε, i.e. Q(u) = Ψ 0 (u) + εΨ 1 (u) + O(ε 2 ), we obtain, from (47), the following two-term equations: To determine how Ψ 1 (u) relates to Ψ 0 (u), we will assume that both c * and Ψ 0 (u) are known, noting that by its construction, Ψ 0 (0) = 0 and Ψ 0 (0) = −c * . Furthermore, from (49), we have that thus, by multiplying by the integrating factor we obtain Noting that Λ(u) ∼ c * as u → 0 + , we determine that as c → c * − ,

Comparison between asymptotic approximations and numerical solutions
With the asymptotic approximations of the heteroclinic orbit Q(u) determined when κ 1, we now compare numerically-obtained solutions of (8) for specific choices of reaction kinetics and diffusivities. Another common generalisation of the Fisher-Kolmogorov model is the Newman reaction-diffusion (NRD) model [4,26], in which The NRD model has an explicitly solvable travelling wave front for a finite critical wavespeed c * , which can be written using the heteroclinic trajectory From (51), we have that implying, from (53), that as c → c * − , where B(p, q) and B x (p, q) are the complete and incomplete Beta functions, respectively [27]. By evaluating (54) at u = 0, we obtain the approximate relationship between κ and c as c → c * − : While (57) cannot be explicitly solved for u(z) for all values of σ, we can expand (57) for u → 0 + to determine the shape of the leading edge of the travelling wave: This expansion confirms that u ∝ (−z) 1/(σ+1) for all c < c * . Furthermore, by expanding (57) for u → 1 − , we obtain the far-field behaviour of u(z): In figure 2, we see that these leading-order approximations for Q(u) and κ(c) are in good agreement with the numerical solution of the NRD model for c → c * − . While approximations of u(z) are less accurate near u = 0, we note that (59) is only a leading-order approximation and further agreement can be obtained with the incorporation of higher-order terms in (57). Nevertheless, the far-field approximation as u → 1 − agrees well with the numerical solution of the NRD model. Our final reaction-diffusion model we will consider is the linear diffusion, oscillatory reaction (LDOR) model, in which Like the FDSR model in section 3.4, we do not attempt to provide a physical justification for the LDOR model, but instead note that, since R(0)D(0) = 1, the critical wavespeed of the LDOR model is c * = ∞. We can therefore use (44) to determine the approximate heteroclinic trajectory when c 1: Additionally, we determine that for c 2 (1 − u) 1, i.e. for u far away from 1, the travelling wave has the approximation where kπ is chosen so that u(z) is continuous on [0, 1]. For c 2 (1 − u) 1, u(z) is described predominantly by the boundary layer of the heteroclinic orbit, implying that Finally, since Q(0) ∼ −c −1 , we determine that κ ∼ c 2 , c 1.
As seen in figure 2, these leading-order approximations for u(z), Q(u) and κ(c) are in good agreement with the numerical solution of the LDOR model for c 1. Specifically, we note that (63) is valid for the majority of the u(z) travelling wave, since the boundary layer near u = 1 is O s (c −2 ) in height (see inset of figure 2(b)). Therefore, we conclude that our leadingorder asymptotic approximations for large-κ travelling waves in the RDS model framework are valid for both c * finite and infinite.

Stability of travelling waves
Now that we have determined various travelling wave solutions of the RDS model, as well as approximate solutions for κ 1 and κ 1, it is natural to ask whether these travelling wave solutions are stable. In general, formal stability analyses of travelling waves arising from reaction-diffusion equations (e.g. [28,29]) have only been considered for travelling waves evolving over R. Furthermore, formal proofs concerning travelling waves in RDS models (e.g. [8]) have been limited to existence and uniqueness of specific choices of reaction and diffusion terms. Due to this gap in the theory of stability analysis, as well as being motivated by simulations of real-life applications, we do not attempt to provide a formal stability analysis of these As shown in figure 3, the GPF and NRD models both agree well with their associated travelling wave solutions, suggesting that these travelling waves are stable. However, the FDSR and LDOR models, whose travelling wave solutions evolve over a finite domain in z, do not agree with numerical solutions of the RDS PDE model for all x. This is expected, since both reaction terms in these models do not equal zero when u = 1 and therefore u continues to grow without bound as t increases. However, we do observe that when κ is large, the section of the numerical solution of the PDE where u ∈ [0, 1] agrees well with the travelling wave profile (figures 3(b) and (d)). Therefore, these comparisons of numerical solutions suggest that RDS models with R(1) = 0 produce unstable travelling wave fronts, while RDS models with R(1) = 0 have stable travelling wave fronts. However, we leave a formal proof of these conjectures, and if there are any additional conditions for stability, for future consideration.

Conclusions
In this work, we consider semi-infinite travelling waves arising from a general reaction-diffusion model coupled with a Stefan-like boundary condition at the leading edge of the front. This RDS model employs a nonlinear diffusivity, D(u), as well as a nonlinear, nonnegative reaction term, R(u). Using travelling wave coordinates, we transform the model into a single nonlinear differential equation, whose solution is the heteroclinic trajectory of the travelling wave in the phase plane. Unlike other reaction-diffusion models, semi-infinite travelling waves arise in the RDS model for a range of wavespeeds, as opposed to a single, critical wavespeed. In the limiting regime where the wavespeed is small, i.e. c 1, we obtain a good approximation of the heteroclinic trajectory of the travelling wave, which can also be used to relate the speed of the front to the Stefan-like parameter, κ. We also determine the approximate form of the travelling wave near its leading edge (u → 0 + ), as well as when u → 1 − . The key result of this analysis is that in order for travelling waves in the RDS model to evolve over a semi-infinite domain, both R(u)D(u) and D(u) must be o(u −1 ) for u → 0 + . Furthermore, travelling waves that evolve over a finite interval require that We also perform a similar asymptotic analysis of travelling wave solutions of the RDS model for wavespeeds approaching a critical wavespeed, c * . This threshold wavespeed provides a bound for when semi-infinite travelling waves cease to exist. Based on the behaviour of R(u)D(u) as u → 0 + , we provide the necessary conditions for which c * is finite. For both c * = ∞ and c * finite, we obtain approximations of the heteroclinic trajectory and the relationship between c and κ as c → c * − . To validate these asymptotic approximations, we examine various choices of R(u) and D(u) in the RDS model that are commonly employed in the literature. In all cases, we find that our asymptotic approximations for the heteroclinic trajectory, the relationship between wavespeed and Stefan-like parameter, and the travelling wave itself all agree well with their corresponding numerical solutions. Finally, a comparison between these travelling wave solutions and numerical solutions of the RDS PDE model is made, as a proxy metric for determining the stability of the aforementioned travelling wave solutions. We find that when R(1) = 0, travelling wave solutions agree well with numerical solutions of the RDS PDE model, suggesting stability. Conversely, if R(1) = 0, reaction terms cause u to increase past u = 1 and do not establish stable travelling wave fronts. Nevertheless, for κ 1, there is agreement between the travelling wave solution and the section of the PDE numerical solutions when u ∈ [0, 1].
Further extensions of the RDS model can also be considered. For instance, one could relax the assumption that R(u) 0 for u ∈ [0, 1], as is done with the Allee-type reaction kinetics [29,30]. For this class of reaction kinetics, the asymptotic analysis described in this work is no longer valid, as the heteroclinic trajectory in the phase plane is no longer strictly negative. While other methods have been proposed to mitigate these issues in related reaction-diffusion models [29], it is unclear how these techniques can be used in the RDS modelling framework. In addition to negative reaction kinetics, other extensions can be incorporated into the RDS modelling framework, including nonlocal reactions [31] and multiple-species reaction-diffusion equations [32,33]. Finally, it should be noted that a formal stability analysis of the travelling wave solutions presented in this work (cf. [28,29]) has yet to be performed. We leave these extensions for future consideration.

Acknowledgments
My sincere gratitude and thanks go to Andrew L Krause (University of Oxford) and Robert A Van Gorder (University of Otago) for their feedback and discussions relating to this work.

Appendix A. Numerical solution for the reaction-diffusion Stefan model
To numerically compute solutions of (1)-(3), we extend numerical methods stated in [11] and first approximate the semi-infinite domain (−∞, L(t)] as the finite domain [0, L(t)], provided that L 0 is sufficiently large. Additionally, we transform (1) and numerically solve the PDE using φ. Consequently, multiplying (1)-(3) by D(u(φ)) and using the transformations stated above yields As we assume that L 0 is sufficiently large, we can replace the far-field boundary condition in u with the Dirichlet boundary condition (A.3) that corresponds to u = 1 when ρ = 0.
To compute the numerical solution of (A.2)-(A.4), we must specify values for L 0 , κ and φ(ρ, 0). We obtain numerical solutions of (A.2) on a uniformly-spaced mesh of ρ ∈ [0, 1], i.e. ρ i = iΔρ, i = 0, . . . , N, where Δρ = 1/N. We denote φ(ρ i , t j ) = φ n i and L(t j ) = L j for convenience, where n 1 is the nth Picard iteration estimate at time t j . Therefore, to determine φ i , we use Here, φ p i is the solution of φ i at the previous timestep, t j−1 , and Δt is the timestep. We identify the system (A.5) and (A.6) as a tridiagonal matrix in φ n i at time t j , which we can solve efficiently using the Thomas algorithm. This solution is stored as Φ n ; if max |Φ n − Φ n−1 | < δ, where δ is some user-specified tolerance, then the Picard loop terminates and we proceed to updating the moving boundary for the next timestep. Otherwise, n = n + 1, the solution Φ n is stored as φ n−1 i , and the Picard iteration loop is performed again. From the fixed-boundary PDE, the Stefan-like condition at ρ = 1 is We approximate φ(ρ, t) as φ(ρ, t j ) during the small interval t ∈ [t j , t j+1 ], allowing us to explicitly solve (A.7) to give a closed form approximation for L(t): Therefore, evaluating this expression at t = t j+1 and using a second-order finite difference approximation of ∂ ρ φ(ρ, t j )| ρ=1 , we obtain the approximation With these updated values of L j+1 and Φ n = φ p i , we update t = t + Δt and j = j + 1; we then repeat the computation to integrate through the next time increment. The algorithm terminates when t + Δt > t f , where t f is the user-specified final time. Once sufficient time has passed that the solution settles towards a travelling wave, we expect that dL(t)/dt = c, so we fit a straight line to our numerical estimate of L(t) and use the slope of that line to provide an estimate of c.