A Nonoscillatory Second-Order Time-Stepping Procedure for Reaction-Diffusion Equations

After a theory of morphogenesis in chemical cells was introduced in the 1950s, much attention had been devoted to the numerical solution of reaction-diﬀusion (RD) partial diﬀerential equations (PDEs). The Crank–Nicolson (CN) method has been a common second-order time-stepping procedure. However, the CN method may introduce spurious oscillations for nonsmooth data unless the time step size is suﬃciently small. This article studies a nonoscillatory second-order time-stepping procedure for RD equations, called a variable-θmethod , as a perturbation of the CN method. In each time level, the new method detects points of potential oscillations to implicitly resolve the solution there. The proposed time-stepping procedure is nonoscillatory and of a second-order temporal accuracy. Various examples are given to show eﬀectiveness of the method. The article also performs a sensitivity analysis for the numerical solution of biological pattern forming models to conclude that the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one. As the spatial resolution becomes higher for an improved accuracy, the CN method may produce spurious oscillations, while the proposed method results in stable solutions.


Introduction
As molecular imaging and single cell analysis is advancing our understanding of spatial processes shaping the cellular dynamics, new models of nonlinear dynamics are necessary. Originating in study of organism development, spatial pattern formation has received a large amount of research over the past decade. Among the most studied, the reactiondiffusion (RD) systems are generating patterns that have been shown to represent well morphogenesis. A theory of morphogenesis based on a RD model was initially proposed by Turing [1]. Gierer and Meinhardt [2] were the first to explore pattern formation in biological systems using the RD equation. Subsequently, several equations of RD type have been studied to understand patterning in developmental biology. Some were derived from phenomenological models (Gierer-Meinhardt) while other modeled simple reaction schemes (Schnackenberg trimolecular autocatalytic reactions model [3], Gray-Scott model [4], Brusselator model [5], chlorite-iodide-malonic acid, CIMA model [6]).
where u � [u 1 , u 2 ] T , D � diag[D 1 , D 2 ] T is the diffusion tensor, Δ denotes the Laplace operator, z/z] is the outward normal derivative on the boundary Γ, and f(u) is the reaction kinetics of the system given as After Turing proposed a theory of morphogenesis in chemical cells in 1952 [1], much attention has been devoted to the numerical solution of RD problems of form (1); see [11][12][13][14] and references therein. Most of the numerical methods studied employed finite difference or finite element approximations for the spatial discretization, while some researchers use finite volume and collocation methods. Once the nonlinear reaction terms are treated (linearized or extrapolated), the Crank-Nicolson (CN) method can be applied as a second-order time-stepping procedure. Timestepping procedures are required at each time step to solve a system of linear algebraic equations, which, although sparse, is compute intensive for multidimensional problems. In order to enhance efficiency of time-stepping procedures, one can adopt the alternating direction implicit (ADI) method as in [11][12][13]15]. In particular, Fernandes et al. [12] introduce an ADI extrapolated CN orthogonal spline collocation method for RD problems.
ADI was invented as a perturbation of the CN method by Douglas, Peaceman, and Rachford in 1955 [16][17][18] and has been employed effectively for the calculation of numerical solution of various time-dependent multidimensional problems, either parabolic or hyperbolic [19,20]. ADI reduces a multidimensional problem to multiple easy-to-solve one-dimensional problems, for an extra cost of a splitting error in O(Δt 2 ), where Δt is the time step size. However, the splitting error can be much larger than the sum of spatial and temporal discretization errors, unless the time step size is sufficiently small [21].
On the contrary, the CN method applied for nonsmooth data may introduce spurious oscillations to the numerical solution unless the time step size is sufficiently small to satisfy the maximum principle 12, which has been recognized in the original paper as well [22]. For this reason, whenever a larger time step or a higher spatial resolution is desirable/necessary, the (less accurate, first order) implicit method which is immune to oscillations has been used at least for several initial time steps with nonsmooth initial data [23]. e CN method and its perturbations (such as ADI) must be applied with care when the solution involves fast transitions or sharp edges; in particular, the time step size should be set very small, e.g., Δt � O(Δx 2 ), where Δx is the spatial grid size. In order to overcome the oscillation problem of the CN method applied for linear parabolic problems of nonsmooth data, the authors recently suggested a variable-θmethod in which the time-stepping parameter of the conventional θ-method, θ ∈ [0, 1], was determined based on local oscillatory characteristics of the solution and the data [24].
In this article, we apply the variable-θ method for the numerical solution of two-component nonlinear RD equations, as given in (1). e variable-θ method is a perturbation of the CN method which evolves the solution implicitly at points where the solution shows a certain portent of oscillations and maintains as a similar accuracy as the CN method with smooth data. e proposed method would be an adequate choice of time-stepping procedure for the numerical solution of RD partial differential equations (PDEs) when a larger time step or a higher spatial mesh resolution is desirable. We have performed a sensitivity analysis for the numerical solution of biological pattern forming models to the spatial and temporal grid sizes. It has been observed from various examples that accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one. When the spatial mesh resolution is set high for a higher accuracy, the method allows to keep the temporal resolution moderate or low. e suggested variable-θ method can result in a smooth/stable numerical solution by suppressing possible oscillations, unlike the CN method. e article is organized as follows. Section 2 includes a brief review on the CN method and its spurious oscillatory behaviors, as preliminaries. In Section 3, a variable-θ method is presented for the numerical solution of two-component nonlinear RD equations. We adopt the successive overrelaxation (SOR) method to solve the resulting linear systems at each time level. Section 4 considers a heuristic technique for the choice of the optimal relaxation parameter for SOR. Section 5 gives numerical examples that show effectiveness of the variable-θ method applied to RD problems in 1D and 2D spaces. In Section 6, we summarize our experiments and present conclusions.

Preliminaries
In this section, we present a brief review of time-stepping procedures, for the numerical solution of linear parabolic equations of the form where D > 0 is a diffusion coefficient and f is a reaction/ source term. We also consider difficulties arising when dealing with nonsmooth data (initial values, boundary conditions, and/or the source term).

e θ-Method: Difference Equation.
Let Ω be a rectangular domain in R 2 : Ω � (a x , b x ) × (a y , b y ). By partitioning Ω × J, we obtain the space-time grid points: 2 Complexity x ij , t n ≔ x i , y j , t n , i � 0, 1, . . . , n x , j � 0, 1, . . . , n y , n � 0, 1, . . . , n t , (4) where n x , n y , and n t are prescribed positive integers and Define the discrete domain, the set of the spatial grid points, by and denote the set of boundary grid points by Γ d � Ω d ∩ Γ and the set of interior grid points by Ω 0 where the FD operators are defined as For the temporal derivative zu/zt, a convenient FD approximation can give Expressing the spatial derivative by a weighted average θ ∈ [0, 1] of previous and current time values, we can formulate the θ-method for (3) as where f n+θ ij is either f(x ij , t n+θ ) or θf n+1 ij + (1 − θ)f n ij . A simple algebraic rearrangement of (10) in a vector form becomes where u n � [u n ij ] 0≤i≤ n x , 0≤j≤n y and f n+θ � [f n+θ ij ] 0≤i≤n x , 0≤j≤n y , considered as column vectors. Popular choices of θ ∈ [0, 1] are 0, 1, and 1/2, which are, respectively, the explicit method (the forward Euler method), the implicit method (the backward Euler method), and the semi-implicit method (the Crank-Nicolson method).
(i) Forward Euler method: when θ � 0, algorithm (11) is stable when it satisfies Although the explicit method is efficient for each time step, its stability condition in (12) enforces, choosing a small time step size Δt; it may become less efficient compared with other implicit methods. It is elementary in numerical analysis that when θ ≥ 1/2, the θ-method (10) is unconditionally stable.
(ii) Crank-Nicolson method: when θ � 1/2, (11) can be rewritten as e CN method has been the most popular time-stepping procedure for the numerical solution of parabolic problems because it is stable and of second-order accuracy in both spatial and temporal directions. However, the CN method applied for nonsmooth data may introduce spurious oscillations to the numerical solution unless the algorithm parameters satisfy the maximum principle [22,23]. As analyzed by the authors [24], the undesired oscillations are due to instability involved in the explicit half step of the CN method, the first term in the right side of (13). e variable θ-method proposed in [24] suppresses spurious oscillations, by evolving the solution implicitly (θ ij � 1) at points x ij where the solution shows a certain portent of oscillations or reduced smoothness, and maintains as a similar accuracy as the CN method with smooth data.
(iii) Backward Euler method: when θ � 1, algorithm (11) reads Although the implicit method shows a first-order accuracy in the temporal direction, it never introduces spurious oscillations to its numerical solutions.

Numerical Oscillations of the CN Method.
Although the CN method is unconditionally stable and of second-order Complexity 3 accuracy in both spatial and temporal directions, it may introduce spurious oscillations into the numerical solution for nonsmooth data. For simplicity, consider a homogeneous diffusion equation with discontinuous initial values defined in the one-dimensional (1D) space: of which the exact solution is given by (16) for (x, t) ∈ (0, 1) × [0, T]. Figure 1 depicts the exact and numerical solutions evolved by the CN and the variable-θ method [24], while Figure 2 compares the numerical solutions at T � 1.0 obtained by the three numerical methods. e numerical solutions are obtained with Δt � 0.01 and Δx � Δy � 0.025. As one can see, spurious oscillations occur along with the step discontinuities in the numerical solution of the CN method, as shown in Figure 1(b). On the contrary, as in Figures 1(c) and 2, the variable-θ method results in an accurate numerical solution without any oscillations. One should notice from Figures 1(a) and 2 that the implicit method is immune to spurious oscillations, but its error is considerably large due to a first-order discretization error in the temporal direction. Although the CN method reveals spurious oscillations, it is quite accurate far from discontinuities. e variable-θ method results in numerical solutions which are smooth as for the implicit method and accurate as for the CN method associated with smooth data. e variable-θ method is a hybrid time-stepping procedure that is based on the CN method (θ � 1/2) and alternately using the implicit method (θ � 1) at points, where the numerical solution shows a certain portent of oscillations or reduced smoothness (the wobble set).

A Variable-θ Method for Two-Component Nonlinear Problems
is section introduces an effective time-stepping procedure for the numerical solution of two-component RD problems (1).

Linearization through Extrapolation.
Once the spatial derivatives are approximated by second-order finite difference schemes, as in Section 2.1, the semidiscrete problem for (1) is formulated as zu zt Let numerical solutions be obtained up to the nth time level, n > 0. For the numerical solution in the (n + 1)th level, we first extrapolate numerical solutions in the two previous levels to approximate the solution at the midpoint t n+1/2 : See [12], for details of second-order extrapolation schemes for n ≥ 0. en, the θ-method for the two-component RD problem reads: e linearized problem (19) can be resolved by solving for two separate components u n+1 � [u n+1 1 , u n+1 2 ] T . Each component in (19) can be formulated as follows: where u, D, and θ denote, respectively, u k , D k , and θ k , for k � 1 or 2, and f n+1/2 is a known source term. e θ-method (20) can be rewritten in a vector form as We present here the main steps of variable-θ method for a nonoscillatory solution of (20); a complete study of the method for diffusion equation was published in [24].

e Variable-θ Method.
e method begins with defining the wobble set, the set of wobble points, as a collection of the grid points where the solution has high fluctuations so that the implicit method (θ � 1) should be applied for the numerical solution not to develop oscillations.
One can easily verify that numerical oscillations of the CN method occur when its explicit half step produces spurious oscillations. Such nonphysical oscillations may happen particularly when the time step size Δt is larger than the stability limit of the explicit scheme. us, the wobble set may be formed to include points where the explicit half step of the CN method introduces undesired local extrema. It follows from (21) that the explicit half step of the CN method (θ � 1/2) reads

Complexity
Let x ij be an interior grid point and consider the four partial directions (made with eight vicinal points of x ij ): four directions having 0°, 45°, 90°, and 135°from the positive x-direction. When spurious oscillations are observed in at least one direction, we select the point x ij as a wobble point.
Define an index function for local extrema (idxt) as Let P, Q, and R be point indices and define iswb(P, Q, R, n) � 1, if idxt u n, * P , u n, * Q , u n, * R ≠ 0 and idxt u n, * P , u n, * Q , u n, * R + idxt u n P , u n Q , u n R < 4, 0, otherwise.
en, the wobble set (for the computation of e function iswb selects candidates for the wobble set from local extrema satisfying idxt(u n, * P , u n, * Q , u n, * R ) ≠ 0; however, the condition idxt u n, * P , u n, * Q , u n, * R + idxt u n P , u n Q , u n R < 4, excludes cases where a strict extremum in u n becomes a strict extremum in the same sense in u n, * . us, the wobble set (25) is the set of interior grid points x ij where u n, * ij becomes a local extremum while u n ij is either a nonextreme or an extreme in the opposite sense, for at least one of four partial directions.
Having the wobble set, the parameter θ for the computation of u n+1 can be assigned pointwise: us, the variable-θmethod for (20) can be formulated as or, in a vector form after grouping variables: e variable-θ method is analyzed for its numerical stability and accuracy and verified for various examples [24]. It results in nonoscillatory numerical solutions of which the accuracy is almost second-order in time.

Remark 2.
e ADI procedure was also applied to (19) for which the initial values show sharp transitions. It has been Complexity 5 observed that ADI may introduce undesirable/discontinuous features to its solution unless the time step size is sufficiently small, i.e., Δt � O(Δx 2 , Δy 2 ). e main problem with ADI is that the diffusion becomes anisotropic, i.e., faster in the coordinate directions. When ADI is applied to the variable-θ formulation of (19), the anisotropic features are reduced significantly. However, it requires to set Δt small enough for a reliable numerical solution. e algebraic system in (29) will be solved by applying the SOR method, with its initial value at the time level t n+1 being set as In particular, SOR converges quite fast for an appropriate choice of the relaxation parameter ω.
In the following, we will consider how to tune the optimal relaxation parameter ω for SOR.

The Optimal SOR Parameter ω
In this section, we will try to find a relaxation parameter which is heuristically optimal. Let us begin with the 2D algebraic system of (21) with θ � 1/2: where L � I + (Δt/2)DA, r � (I − (Δt/2)DA)u+ Δtf n+1/2 , and m > 0 is the dimension of the algebraic system. It is known that the optimal relaxation parameter for the SOR method can be determined as ( [25], Section 4.3) where ρ(T J ) is the spectral radius of the Jacobi iteration matrix T J . For simplicity, assume that the problem is defined on the unit square with a Dirichlet boundary condition. We further assume that the domain is partitioned into N × N grids so that h � Δx � Δy � 1/N. en, the eigenvalues of the second-order 5-point FD coefficient matrix A read ( [25], Section 6.5) and therefore the eigenvalues of L can be formulated as for 1 ≤ k, ℓ ≤ N − 1. Note that the diagonal element of L is So, the eigenvalues of the Jacobi iteration matrix T J are given as ) .

(36)
In order to find the maximum of |λ k,ℓ (T J )|, we first obtain for some c 1 > 0. Here, we have used h � 1/N and the approximation cos(x) ≈ 1 − x 2 /2. Now the spectral radius of T J reads Assuming that c 1 h 2 < h 2 /(2 DΔ t) < 1, we finally obtain for some c 2 > 0. It follows from (32) and (39) that the optimal SOR parameter ω Δt,h , corresponding to the spatial grid size h and the time step size Δt, can be written as for some c 0 > 0. e constant c 0 can be found experimentally from a selected set of (Δt, h), as summarized in the following: (a)Determine ω Δt 0 ,h 0 for prescribed grid sizes Δt 0 , h 0 , heuristically.

Numerical Experiments
In this section, we present numerical experiments which show effectiveness of the variable-θ method. e algorithms are implemented in MATLAB and carried out on a desktop computer of Intel Xeon CPU E5-1620 3.60 GHz processor.
To solve the algebraic system at each time level, the SOR method is employed with the near-optimal parameter ω calculated as in (40), with c 0 being estimated with a small grid problem. e SOR iteration is stopped when the maximum difference of consecutive iterates becomes smaller than a tolerance ε � 10 − 6 : 6 Complexity e L ∞ -error E ∞ [t n ], measured at t � t n , is defined as follows: where u is the exact solution.
5.1. One-Component RD System. To investigate accuracy of the variable-θ method, we consider the diffusion problem (15) studied earlier in Section 2.2. For a comparison purpose, we have implemented not only the θ-methods and the variable-θ method but also the implicit predictor-corrector (0, 2)-Padé (IPC-[0,2]) method [26] and Lawson and Morris (LM) local extrapolation method [23]. Table 1 presents the L ∞ -error E ∞ [T] at T � 1.0 when the five methods are applied for the numerical solution of (15), with various Δt and Δx. As one can see from the table, the CN method resolves its numerical solution poorly (due to oscillations), except for the case the method satisfies the maximum principle. On the contrary, the variable-θ method results in a second-order accuracy, with its errors being smallest among all the methods for most cases. e new method is a hybrid time-stepping procedure which assembles merits from the CN method (high-accuracy) and the implicit method (smoothness). Now, consider a nonlinear RD problem of the form with the boundary and initial values, as given in (15). Figure 3 presents the numerical solutions evolved by the implicit method, the CN method, and the variable-θ method, when Δt � 0.01 and Δx � 0.025 (the mesh is the same as the one selected in Figure 1). Similar to the linear problem in Figure 1, spurious oscillations are introduced into the numerical solution of the nonlinear problem by the CN method only. It should be noticed that spurious oscillations of the CN method are damped out much faster for the nonlinear problem than the linear problem, which is due to the reaction kinetic term f(u) � u(1 − u). For the nonlinear problem, It seems that the oscillations at early time steps do not affect the solution at later steps much. is observation explains a partial reason that the second-order CN method has been popular for the numerical solution of PDEs in mathematical biology. However, for other applications, such spurious oscillations at early moments may alter the numerical solution significantly so as for the CN method to be unstable; see Figure 4 below. It is important to develop an effective algorithm which can suppress spurious oscillations for convenient choices of algorithm parameters; the variable-θ method is effective and stable.

Two-Component Nonlinear RD Systems.
Two-component RD systems enable to explain a much wider range of phenomena than their one-component counterparts. Many two-component models have been developed and numerically verified for dynamical patterning behaviors in biology and chemistry. In this section, we consider two-component models interested in the literature of biology and chemistry, to verify effectiveness of the variable-θ method.

Gray-Scott Model in 1D.
We apply the numerical methods for the numerical solution of the Gray-Scott model [4,27] defined as (1) associated with the following reaction kinetics: for any constants F and k. Let Ω � (0, 1). We assign two sets of model constants and initial and boundary conditions as follows [28]: where (46b) is a mid-pulse initial condition, and D � 10 − 4 , 5 · 10 − 5 T , F � 0.025, k � 0.0544, (47a) x ∈ (0, 1), where (47b) is a left-pulse initial condition. In Figure 5, we present the propagation of the numerical solution of u 2 by the variable-θ method associated with (46a)-(46c) at the grid sizes Δt � 0.01 and Δx � 0.004 over 0 ≤ t ≤ T � 2000. e initial midpulse splits in early moments to travel in both directions (self-replication of the pulse). As each of the pulses travels, it becomes thicker (bigger) up to a certain width and begins to replicate itself recursively. Complexity Figure 6 depicts the propagation of the numerical solution of u 2 by the variable-θ method associated with (47a)-(47c) over 0 ≤ t ≤ T � 5000, at the same resolution, as in Figure 5. One can clearly observe a traveling pulse which begins from the left edge point and reflects whenever it hits the boundary, due to the no-flux boundary condition (47c).
To investigate bifurcation in the RD system and our method numerical accuracy, we present numerical solutions of the wave-splitting problem (46a)-(46c) obtained with various spatial and temporal grid sizes, as shown in Figure 7. e image I kℓ represents the numerical solution obtained with the mesh resolution (Δt, Δx) � (10 − k , 0.01/2 ℓ− 1 ). For example, the image I 23 is associated with the mesh resolution (Δt, Δx) � (1/100, 1/400). One can easily point out from the images that the spatial resolution alters the numerical solution dramatically even with halved spatial grid sizes (compare the images horizontally), while the temporal resolution affects little the numerical solution even with oneorder smaller temporal step sizes (compare them vertically). e main reason for such a sensitivity to the spatial resolution is that the RD does not have as much time before growing to reach the margins of the mesh in low spatial resolutions (of large Δx's). When this happens, the RD pattern typically deteriorates and it does not travel in an appropriate speed nor reaches a condition to replicate itself on time, see ( [29], Section 4.2) for similar observations. We summarize the experiments with the Gray-Scott model in 1D as follows: (i) Accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one. (ii) us, it is crucial to set a high spatial resolution (small Δx's) for a desirable accuracy.
(iii) As the spatial resolution becomes higher, the CN method may more likely produce spurious oscillations, while the variable-θ method results in stable solutions.

Gierer-Meinhardt Model.
e Gierer-Meinhardt model [2] is (1) defined in Ω � (− 1, 1) × (− 1, 1) ⊂ R 2 with the following reaction kinetics and parameters: for which various numerical methods have been developed [12,31,32]. We cast the experiment employing coefficients and the initial condition used in [32]: e initial values are depicted in Figure 9. In this section, we restrict our attention to the dynamics of u 1 of the model. In order to investigate effectiveness of the variable-θ method and oscillatory behaviors of the CN method as well, we have carried out numerical experiments for the Gierer-Meinhardt model with a relatively low spatial resolution. Figure 4 presents numerical solutions at two different times for u 1 of the Gierer-Meinhardt model with the spatial resolution Δx � Δy � 1/32. When the time step size is set Δt � 0.05, the variable-θ method evolves the numerical solution as shown in Figures 4(a) and 4(d), for which the final steady-state pattern is the same as that in [32]. On the contrary, with the same Δt � 0.05, the CN method has produced a quite different pattern as in Figures 4(b) and 4(e), due to the nonsmooth initial values (Figure 9). However, the CN method can recover the correct steady-state pattern when it runs with Δt � 0.005, as depicted in Figures 4(c) and 4(f ). For a similar accuracy, the variable-θ method (taking 170 s) is about 7 times more efficient than the CN method (taking 1242 s).
We summarize our experiments with the Gray-Scott and Gierer-Meinhardt models in 2D as follows: (i) e variable-θ method shows the same accuracy as the CN method for problems of smooth data (ii) For nonsmooth data, the variable-θ method evolves a smooth solution for all choices of Δt, while the CN method introduces spurious oscillations to alter the solution unless the time step size is sufficiently small (iii) When a large time step size is desirable, the suggested method is a few times more efficient than the CN method for a similar accuracy    2 (x, y, 0).

Remark 3.
Although the variable-θ method can employ larger time step sizes than the CN method to get stable numerical solutions for problems of nonsmooth data, one may not set the time step size too large, due to an accuracy issue rather than the stability issue. Furthermore, for nonlinear problems, the overall stability of the numerical algorithm can be determined by not

Complexity 13
only grid sizes but also numerical schemes including methods of dealing with the nonlinear terms. Figure 10 presents numerical solutions and their errors of u 1 for the Gray-Scott model (48)-(50) at T � 1.0 by the variable-θ method, varying Δx � Δy with fixed Δt � 0.25. Compared with Figure 8, the solutions show stability and a good accuracy, although the time step size is as large as Δt � 0.25. As shown in the bottom line in Figure 10, all three spatially different cases show the same level of errors since the entire errors are dominated by temporal direction errors. We conclude from this example that grid sizes in both temporal and spatial directions would not significantly affect the stability of the proposed method when the initial condition is smooth and the nonlinearity is not severe.
As an example of nonsmooth data and severe nonlinearity, we select the Gierer-Meinhardt model (51) and (52) to simulate with large temporal step sizes. When Δt ≥ 0.2, the proposed algorithm introduced a rapid decay of solution values independently of the spatial grid size, so that the pattern is not formed appropriately. We believe that it is due to the error incorporated with the reaction term (19) when u n+1/2 is approximated by the extrapolation scheme (18). However, when Δt ≤ 0.1, our method produces stable solutions for all choices of spatial grid sizes. Figure 11 presents numerical solutions of u 1 for the Gierer-Meinhardt model at T � 340 by the variable-θ method with fixed Δt � 0.1 and various Δx � Δy. Note that for Gierer-Meinhardt model, the pattern forming is slow down as the spatial grid size becomes smaller, as shown in (I 2 ) and (I 3 ) of Figure 11; this tendency has been observed for all other choices of Δt ≤ 0.1.
is is another example that accuracy of the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.

Conclusions
e Crank-Nicolson (CN) method has been a popular second-order time-stepping procedure for the numerical solution of systems of nonlinear RD equations. However, the CN method may introduce spurious oscillations for nonsmooth data unless the time step size is sufficiently small. We have studied a nonoscillatory time-stepping procedure for RD equations, called a variable-θ method, as a perturbation of the CN method. In each time level, the new method detects points of potential oscillations and resolves the solution applying the implicit method locally at those points. e proposed timestepping procedure has proven nonoscillatory and having a second-order temporal accuracy, although the initial conditions are nonsmooth. Various examples have been considered to show effectiveness of the method. We also have performed a sensitivity analysis for the numerical solution of biological pattern forming models to conclude that the numerical solution is much more sensitive to the spatial mesh resolution than the temporal one.
Data Availability e experiment data and figures used to support the findings of this study are included within the article, and the MATLAB codes for the experiments are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.