A Mass- and Energy-Conserving Numerical Model for a Fractional Gross–Pitaevskii System in Multiple Dimensions

: This manuscript studies a double fractional extended p -dimensional coupled Gross– Pitaevskii-type system. This system consists of two parabolic partial differential equations with equal interaction constants, coupling terms, and spatial derivatives of the Riesz type. Associated with the mathematical model, there are energy and non-negative mass functions which are conserved throughout time. Motivated by this fact, we propose a ﬁnite-difference discretization of the double fractional Gross–Pitaevskii system which inherits the energy and mass conservation properties. As the continuous model, the mass is a non-negative constant and the solutions are bounded under suitable numerical parameter assumptions. We prove rigorously the existence of solutions for any set of initial conditions. As in the continuous system, the discretization has a discrete Hamiltonian associated. The method is implicit, multi-consistent, stable and quadratically convergent. Finally, we implemented the scheme computationally to conﬁrm the validity of the mass and energy conservation properties, obtaining satisfactory results.


Introduction
Fractional calculus emerged as a generalization of conventional calculus. The classical differential and integral operators proposed by Leibniz more than three centuries ago have been extended to the fractional scenario in various different forms. In the way, various applications have been proposed to within mathematics, such as Cauchy problems with Caputo Hadamard fractional derivatives [1], the synchronization for a class of fractionalorder hyperchaotic systems [2], the inequality estimates for the boundedness of multilinear singular and fractional integral operators [3], the analysis of unstructured-mesh Galerkin finite element method for the two-dimensional multi-term time-space fractional Bloch-Torrey equations on irregular convex domains [4], and the design of numerically efficient and conservative model for a Riesz space-fractional Klein-Gordon-Zakharov system [5], among various other interesting problems [6].

(4)
This system is particularly interesting due to the fact that it has quantities which are conserved through time, one of them being the energy function. More precisely, the total energy at time t ∈ [0, T] is provided as Definition 2. For each z ∈ C, z represents the complex conjugate of z. Let L x,2 (Ω) denote the set of all functions f : Ω → C with the property that, for each t ∈ [0, T], the function f (·, t) : Ω → C belongs to L 2 (Ω). For each pair f , g ∈ L x,2 (Ω), the inner product of f and g is the function of t defined by In addition, if u = (u 1 , u 2 , . . . , u p ) and v = (v 1 , v 2 , . . . , v p ) are vector functions such that u i , v i ∈ L x,2 (Ω) for each i ∈ I p , then we set u, v x = u 1 , v 1 x + u 2 , v 2 x + . . . + u p , v p x . Let f x,2 = f , f , for each f ∈ L x,2 (Ω). In similar fashion, we can define the space L x,p (Ω). If f ∈ L x,p (Ω), then let The next result summarizes the property of conservation of energy satisfied by (4). x,4 + β 12 ψ 1 ψ 2 2 x,2 + 2λ Re ψ 1 , ψ 2 x , for each t ∈ (0, T). Moreover, if ψ 1 and ψ 2 satisfy (4), then the function E (t) is a constant.
Proof. The proof of this result was given by Serna-Reyes and Macías-Díaz [20].
The energy density of the system (4) is defined as the function in the integrand of the expression of E (t). More precisely, the energy density of the physical model (4) is given by the function H(ψ 1 , On the other hand, the total mass of the continuous system (4) at the time t is defined as M(t) = M 1 (t) + M 2 (t), for each t ∈ (0, T). Here, we let M 1 (t) = ψ 1 2 x,2 and M 2 (t) = ψ 2 2 x,2 , for each t ∈ (0, T).
Theorem 2 (Mass conservation). If ψ 1 and ψ 2 are solutions of (4), then the function M(t) is non-negative and constant. Moreover, if λ = 0, then the functions M 1 (t) and M 2 (t) are also non-negative and constant.
Proof. Beforehand, recall that the square-root operator of − α is the operator α/2 , for each α ∈ (1,2]. This means that − α u, v x = α/2 u, α/2 v x , for any functions u, v ∈ H 2 (Ω T ). Notice now that the following identities hold: These inner products are real numbers, except the first one which is purely imaginary. On both sides of the first equation of (4), obtain the inner product with ψ 1 and take the imaginary parts. In addition, calculate the inner product on both sides of the second differential equation of (4) with ψ 2 , and take also the imaginary parts on both sides. Then, Notice that the individual masses M 1 (t) and M 2 (t) are constant when the internal atomic Josephson junction is equal to zero, as mentioned in the conclusion of this result. Moreover, if we add these two identities, it readily follows that dM(t) dt = λ Im ψ 1 , ψ 2 x + ψ 1 , ψ 2 x = 2λ Im(Re ψ 1 , ψ 2 x ) = 0, ∀t ∈ (0, T). (18) Conclude that the total mass of the system is conserved throughout time, as desired.
In the present work, we design a finite-difference method to solve the system (4), in such way that the main structural properties of the solutions of that continuous system are also satisfied in the discrete domain. More precisely, we design a numerical method which is able of preserving the discrete energy and mass of the system, and such that the solutions satisfy similar boundedness properties. The scheme is thoroughly analyzed for consistency, stability, and consistency, and we illustrate these features through some computer simulations. As it turns out, the proposed numerical model has an easy implementation, as shown by the code in Appendix A.

Numerical Model
In the following, we propose a linearly-implicit numerical method to approximate the solutions of the fractional system (4), along with the energy and mass functions. To that end, the concept of fractional centered differences is of utmost importance.
Definition 3 (Ortigueira [21]). Let f : R → R be any function, h > 0, and α > −1. The fractional-order centered difference of f of order α at x is defined as where g (α) Lemma 1 (Wang et al. [22]). Let 0 < α ≤ 2 and α = 1. Then, Lemma 2 (Wang et al. [22]). Let 0 < α ≤ 2 and α = 1. If f : R → R has integrable derivatives up to order 5, then, for almost all x, Next, we present the discrete notation used to approximate the solutions of the fractional problem (4). We let h i , τ ∈ R, for each i ∈ I p . Moreover, we assume that N = T/τ and M i = (b i − a i )/h i are natural numbers, for each i ∈ I p , and fix uniform partitions of [a i , b i ] and [0, T] given, respectively, by t n = nτ, ∀n ∈ I N .
Let J = ∏ p i=1 I M i −1 and J = ∏ p i=1 I M i , and let ∂J represent the boundary of the mesh J. Define x j = (x 1,j 1 , . . . , x p,j p ), for each multi-index j = (j 1 , . . . , j p ) ∈ J. The symbol (u n j , v n j ) denotes a numerical estimate to (ψ 1 (x j , t n ), ψ 2 (x j , t n )), for each (j, n) ∈ J × I N .
Definition 4. Define the following, for w = u, v and α ∈ (0, 1) ∪ (1, 2], (j, n) ∈ J × I N−1 : and δ (α) . , x p,j p , t n ). (33) Moreover, for the sake of convenience, we agree that The finite-difference method to approximate the solution of (4) on Ω T is the discrete system with initial-boundary conditions iδ (1) where (j, n) ∈ J × I N−1 . Notice that this scheme is an implicit three-step method. Moreover, the approximations at the time t 0 are provided by the initial data u 0 On the other hand, the discrete model requires knowledge of fictitious approximations at the time t −1 . To avoid this shortcoming, notice that the initial conditions µ (1) j , respectively, for each j ∈ J. Substituting these identities into the recursive difference equations of (36) when n = 0, we obtain the explicit formulas and For the sake of convenience, the stencil of the scheme (36) is provided in Figure 1.
We study the existence and uniqueness of solutions of the finite-difference method (36). To that end, we require some additional nomenclature and technical results. To start with, we let h = (h 1 , . . . , h p ) and h * = ∏ p i=1 h i , and fix the spatial mesh {x j } j∈J ⊆ R p . Let V h be the complex vector space of all grid functions which vanish at the boundary of the mesh. For any w ∈ V h and j ∈ J, convey that w j = w(x j ).

Definition 5.
Define the product ·, · : V h × V h → C and the norm · 1 : V h → R by The Euclidean norm induced by the inner product ·, · is denoted by · 2 . x j+2 x j−2 . . . × × Figure 1. Stencil of the numerical model (36) at the time t n in one spatial dimension. The black circles denote the known solutions at the times t n−1 and t n . Meanwhile, the crosses represent the unknown approximations at the time t n+1 . For the sake of convenience, we agree that x j i = x 1,j i .
In the following, we represent the solutions of the finite-difference method (36) by (u n , v n ) N n=0 , where we agree that u n = (u n j ) j∈J and v n = (v n j ) j∈J , for each n ∈ I N . We use the following lemma to prove that solutions of (36) exist.
Lemma 3 (Browder fixed-point theorem [23]). Let (H, ·, · ) be a finite-dimensional innerproduct space, let · : H → H be the norm induced by the inner product, and suppose that g : H → H is continuous. Assume that there exists β > 0 such that Re g(z), z > 0, for all z ∈ H with z = β. Then, there is z * ∈ H, such that g(z * ) = 0 and z * ≤ β.

Theorem 3 (Solubility).
For any set of initial conditions, the discrete system (36) is solvable.
Proof. The proof requires mathematical induction. Notice that the approximations at the times t 0 and t 1 are provided by the initial profiles and the explicit systems (37) and (38), respectively. Thus, let us suppose that the numerical solutions of the discrete model (36) at the times n − 1 and n are known, for some n ∈ I N−1 . Rearranging the recursive equations of the method, we can readily check that the following equations hold: and for each η ∈ V h , ν ∈ V h , j ∈ J. After some calculations, we readily reach that Using now the Cauchy-Schwarz inequality and bounding from below, it is easy to check that the following hold: If we let β = (u n−1 , v n−1 ) + λτ (v n , u n ) + 1 in the previous lemma, it follows that there exists (u n+1 , v n+1 ) ∈ V h × V h which satisfies the discrete system (36). The conclusion of this result readily follows now by induction.

Structural Properties
The present section is devoted to showing that the discrete model (36) satisfies physical properties which are similar to those satisfied by the continuous system (4). More precisely, we propose numerical energy and mass functional associated to the scheme (36) which are conserved conditionally. Various technical lemmas are required to that end. [24]). For each i ∈ I p and k ∈ I 2 , there exists a unique positive operator

Lemma 5.
The following equations hold, for each n ∈ I N−1 :

Re iδ
Proof. To establish (48), notice that the following identities trivially hold: To prove the identity (49), observe that Re Vµ which establishes the second identity. We prove now the formula in part (51). To that end, notice that Finally, we must point out that the proofs for the identities (50) and (52) are similar to the proofs of (49) and (51), respectively, and, for that reason, we omit them here. Lemma 6. The following equations hold, for each n ∈ I N−2 and k ∈ I 2 : Re( v n , 2δ (1) Proof. Applying distributivity of the inner product and rearranging terms, whence the first identity readily follows. To prove the second identity, we use once again the distributivity of the inner product and Lemma 4 to see that which readily establishes the identity. Finally, using again distributivity and rearranging terms conveniently, Re v n , 2δ (1) t u n + u n , 2δ (1) This last chain of identities readily establishes the last property of the lemma.
The next theorem proves the existence of energy-like invariants for (36).
Theorem 4 (Energy conservation). Let (u n , v n ) N n=0 be a solution of (36), and define for each n ∈ I N−1 . Then, δ t E n = 0, for each n ∈ I N−1 .
Proof. Denote the left-hand sides of the two difference equations in (36) by Θ n j and Φ n j , respectively, for each j ∈ J and n ∈ I N−1 , and define n=0 is a solution of (36), calculating the inner product of Θ n with 2δ (1) t u n and of Φ n with 2δ (1) t v n , using the identities above and collecting terms, we reach In similar fashion, notice that the following hold, for each n ∈ I N−1 : Adding these two last identities, we readily obtain that δ t E n−1 = 0, for each n ∈ I N−1 . The conclusion of this result is readily reached now using mathematical induction.
Motivated by this result, the quantities E n are called the discrete total energy of the system (36) at the time t n , for each n ∈ I N−1 . In addition, we define the discrete energy density of the system (36) as for each (j, n) ∈ J × I N−1 .

Lemma 7.
The following equations hold, for each n ∈ I N−2 and k ∈ I 2 : Im Vµ Im Dµ Im iδ Proof. Identity (a) readily follows from the fact that which obviously yields zero. On the other hand, applying distributivity of the inner product on both factors, we obtain Im iδ which establishes Property (f). Finally, notice that distributivity also yields Im v n , 2µ (1) t u n + u n , 2µ (1) which shows the validity of (g). The proofs of Properties (b)-(e) are similar to that of (a).
For that reason, we omit the proofs of the remaining identities.
Theorem 5 (Mass conservation). If (u n , v n ) N n=0 is a solution of (36), then δ t M n = 0 for each n ∈ I N−2 , where is the discrete total mass of the system at the time t n . The discrete individual discrete masses as the quantities and they are conserved if λ = 0.
Proof. Let Θ n and Φ n be as in the proof of Theorem 4, for each n ∈ I N−2 . Calculating the inner product of Θ n and Φ n with 2µ (1) t u n and 2µ (1) t v n , respectively, using the identities in the previous lemma. Collecting terms, we note that and for each n ∈ I N−2 . Observe that, indeed, the individual discrete masses are conserved if λ = 0. On the other hand, adding the last two identities, we readily check that for each n ∈ I N−2 , which is what we wanted to prove.
After this last result, we refer to the quantities M n as the total mass of the discrete model (36) at the time t n , for each n ∈ I N−1 . Moreover, we define the density of mass at the point (x j , t n ) as the number for each (j, n) ∈ J × I N−1 . It is worth noticing that the quantities E n and M n in the last theorems are clearly the discrete counterparts of the total energy and the total mass of associated to the continuous system (36), respectively.
Corollary 2 (Mass non-negativity). Let (u n , v n ) N n=0 be a solution of (36), and suppose that |λ|τ < 1. Then, the quantities M n are non-negative and constant.
Proof. It only remains to prove that the quantities M n are non-negative. Using Young's inequality, we obtain that for each n ∈ I N−1 . This obviously completes the proof.
The proof of the following result is similar to that of Corollary 2.

Numerical Properties
To start with, we introduce the continuous functional for each (x, t) ∈ Ω T , and the next discrete functions, for each (j, n) ∈ J × I N−1 : and Using this nomenclature, we agree that for each (x, t) ∈ Ω T and (j, n) ∈ J × I N−1 .
Theorem 6 (Consistency). If ψ 1 , ψ 2 ∈ C 5,4 x,t (Ω T ), then there exist constants C, C , C > 0 which are independent of h and τ, such that, for each (j, n) ∈ J × I N−1 , the following hold: Proof. Use Taylor's theorem, Lemma 2 and the conditions of regularity, there are constants C 1,k , C 2,k,i , C 3,k , C 4,l ∈ R which are independent of h and τ, for i ∈ I p , k ∈ I 2 and l ∈ I 3 , such that the next hold, for each (j, n) ∈ J × I N−1 : The first inequality of the conclusion follows from the triangle inequality. The proof of the second inequality is similar. Finally, to prove the last inequality, observe firstly that the regularity conditions on ψ 1 and ψ 2 assure that there exist constants C 5,k , C 6 ∈ R which are independent of h and τ, for k ∈ I 2 , such that The third inequality of this result readily follows now assuming that τ < 1, using the triangle inequality, and letting C be a suitable constant in terms of C 5,1 , C 5,2 and C 6 .
Lemma 8 (Ferreira [25]). Let (ω n ) N n=0 and (ρ n ) N n=0 be arrays of non-negative numbers for which there exists C 0 ≥ 0 with the property that the following hold: If n ∈ I N , then ω n ≤ ρ n e C 0 nτ .
To prove stability, take two solutions of our scheme (36) associated to two sets of initial data. In notation, we assume that (u, v) represents a solution for (36), and that (ũ,ṽ) solves the discrete problem but with initial profiles described by functionsφ 1 ,φ 2 : Ω → C. Concretely, the following discrete initial-boundary-value problem holds: for each (j, n) ∈ J × I N−1 .
Lemma 9. Let = ( n ) N n=0 and ζ = (ζ n ) N n=0 be sequences in V h . If m ∈ I N−1 then τ m ∑ n=1 δ t Im n−1 , ζ n + ζ n−1 , n = Im m , ζ m+1 + ζ m , m+1 − 0 , ζ 1 − ζ 0 , 1 . (103) Proof. After substituting the definition of the operator δ t at the left-hand side of the identity (103) and performing algebraic operations, we obtain that Observe that most of the terms inside the sums cancel out. The conclusion readily follows now. Lemma 10. Let (u, v) and (ũ,ṽ) be solutions of (36) and (102), respectively, and let us define n = u n −ũ n and ζ n = v n −ṽ n , for each n ∈ I N . For each (j, n) ∈ J × I N−1 , define There is a constant C ≥ 0 which is independent of h and τ, such that, for each (j, n) ∈ J × I N−1 , Proof. It is straightforward.

Theorem 7 (Stability).
Assume that (u, v) and (ũ,ṽ) are solutions of the problems (36) and (102), respectively, and define n = u n −ũ n and ζ n = v n −ṽ n , for each n ∈ I N . Assume that the condition 2λ + 4C + 1 τ < 1 (108) is satisfied. Then, there exists a constant C ≥ 0 which is independent of τ and h, such that Proof. Obviously, the next problem is satisfied, for each (j, n) ∈ J × I N−1 : where P n j and Q n j are as in Lemma 10. Then, there is C ≥ 0 independent of h and τ, with the property that, for each n ∈ I N−1 , max{ P n 2 2 , Q n 2 2 } ≤ C n−1 2 Obtain the inner product of the first vector Equation of (110) with 2µ (1) t n and calculate the imaginary parts. Meanwhile, in the second equation, obtain the inner product with 2µ (1) t ζ n and calculate imaginary parts. By Theorem 5, t ζ n + 2λ Im n , µ t ζ n , ∀n ∈ I N−1 .
It is readily checked then that that there exists C ≥ 0 independent of τ and h, such that (114) The conclusion of this result is reached now using Lemma 8.
The following result is now straightforward.

Corollary 4 (Uniqueness).
For any initial data, the problem (36) has a unique solution.
We establish next the convergence properties of our finite-difference scheme. To that end, we assume that (ψ 1 , ψ 2 ) represents a pair of sufficiently smooth solutions of the initial-value problem (4). As a consequence, Obviously, ζ n j and ξ n j represent here the local truncation errors.

Numerical Simulations
The present section is devoted to describe our computational implementation of the finite-difference scheme (36), and to provide numerical simulations that illustrate the capability of the scheme to preserve the energy and the mass of the system. Our description of our implementation and the simulations are carried out in the one-dimensional case, that is, when p = 1. Moreover, for simplification purposes, let a = a 1 , b = b 1 , h = h 1 , and M = M 1 . Throughout this section, we represent the approximations in vector form as For computational purposes, we define next various matrices, all of them with sizes equal to (J + 1) × (J + 1). To start with, we let A n u = I + D − u + 1 2 τG (α 1 ) − τE n u and B n u = D + u − 1 2 τG (α 1 ) + τE n u , for each n ∈ I N , where I, D ± u and E n u are diagonal matrices defined through their diagonal entries by the formulas and E n u = diag 0, β 11 |u n 1 | 2 + β 12 |v n 1 | 2 , β 11 |u n 2 | 2 + β 12 |v n 2 | 2 , . . . , . . . , β 11 |u n J−1 | 2 + β 12 |v n J−1 | 2 , 0 , ∀n ∈ I N , and for each α ∈ (1, 2], we define For each n ∈ I N , we also introduce the matrices A n v = I + D − v + 1 2 τG (α 2 ) − τE n v and B n v = D + n − 1 2 τG (α 2 ) + τE n v , for each n ∈ I N , where the new matrices D ± v and E n v are diagonal, and they are defined as and E n v = diag 0, β 22 |v n 1 | 2 + β 12 |u n 1 | 2 , β 22 |v n 2 | 2 + β 12 |u n 2 | 2 , . . . , . . . , β 22 |v n J−1 | 2 + β 12 |u n J−1 | 2 , 0 , ∀n ∈ I N , Observe that the approximations at the times t 0 of the discrete model are provided exactly through the initial conditions in (36). Moreover, the approximations at the time t 1 are provided explicitly through Formulas (37) and (38). Meanwhile, the recursive equations of the discrete model can be alternatively expressed in terms of the nomenclature introduced in this work as where the (J + 1)-dimensional vectors u n * and v n * are given by u n * = 0, u n 1 , u n 2 , u n 3 , . . . , u n J−1 , 0 , ∀n ∈ I N , It is worth pointing out that our computer implementation in Appendix A is based on the recursive system (134).

Example 1.
In the following, we use the constants in Table 1. Meanwhile, the parameters α 1 and α 2 change values between simulations, and V(x) = 1 2 x 2 , for each x ∈ Ω. On the other hand, for each x ∈ Ω. Throughout, we fix the computer parameters τ = 0.05 and h 1 = 0.1. In our simulations, we consider three cases for the differentiation orders α 1 and α 2 , namely, Case 1. α 1 = 1.1 and α 2 = 1.9.  Under these circumstances, Figures 2-4 provide the approximate solutions for the problem (4), using the values of Cases 1-3, respectively. In each figure, the graphs provide the behavior of the following , versus (x, t): (a) Re ψ 1 ; (b) Re ψ 2 ; (c) Im ψ 1 ; (d) Im ψ 2 ; (e) |ψ 1 |; (f) |ψ 2 |. Figure 5 provides approximations for the energy functionals for these problems. More precisely, the left column of that figure provides approximate graphs of the discrete energy density H versus (x, t), while the graphs of the total energy E versus t are shown on the right column. We used Case 1 (top row), Case 2 (middle row) and Case 3 (bottom row). Obviously, the total energy of the system is approximately constant, as expected from the theoretical results. Finally, Figure 6 provides a similar analysis for the mass functionals of the mathematical model (4). The results show that the mass is approximately constant in all the cases, in agreement with the theoretical results proved in this manuscript.
(139) Before closing this section, it is important to point out that the extension of the Gross-Pitaevskii system to the fractional scenario (4) was proposed as a natural generalization of the classical Bose-Einstein condensates. Within that context and to the best of our knowledge, the theory of fractional operators has not found a solid justification via experimentation. However, the mathematical and numerical investigation of fractional extensions of this system is a difficult problem whose interest circumscribed to mathematical and numerical areas. Nevertheless, the authors hope that some potential applications of the present work will surface in the near future using some well-known applications of fractional derivatives (see Appendix B), and that the mathematical and numerical ideas developed in this work will lead to propose and analyze new computational methods for other physically relevant systems.

Conclusions
In this work, we study a finite-difference method for the coupled system of Gross-Pitaevskii equations with space-fractional derivatives in the Riesz sense. We prove the existence of solutions of our discretization for any set of initial conditions using a complex fixed-point theorem. Inspired by the existence of invariants in the continuous system, we show that some discrete forms of energy and mass functions are constant in time. As its continuous counterpart, we also prove the non-negativity of the mass function and the uniform boundedness of the approximations by imposing suitable numerical con-straints. We present a rigorous numerical analysis for the method, obtaining a second order of consistency in space and time for the solutions. Moreover, we show that the discrete energy and mass densities are consistent approximations of the continuous functionals. We also provide conditions to guarantee the stability and the uniqueness of the numerical solutions, and we prove the convergence of the system. The theoretical results show that the method is quadratically convergent. It is worth pointing out that analyses of stability and convergence were carried out using a discrete version of Gronwall's inequality. A Matlab implementation of the numerical model is provided in the Appendix of this work, and we used to produce approximations to the solutions of our continuous model. The results show that the energy and mass are approximately constant, in agreement with the theoretical results derived in this work, and improving computationally efforts already reported [20,26].
It is worth pointing out that there are already some reports by these same authors on the design of numerical methods to solve double Gross-Pitaevskii-type systems with fractional diffusion, most notably the papers by Serna-Reyes et al. [26] and Serna-Reyes and Macías-Díaz [20]. In both works, the same system of fractional partial differential equations is investigated using a finite-difference methodology, the fractional derivatives are of the Riesz type, and they are discretized using fractional-order centered differences. To start with, the proofs for the existence of solutions of the numerical model are more complicated in the case of the published papers [20,26]. Indeed, the arguments in those cases require using some fixed-point theorem in view of the nonlinear nature of those works. Meanwhile, in the present work, the existence of solutions is readily established using matrix arguments. This is due to the linear nature of the numerical solved. Moreover, the uniqueness in the present case is derived at the same time as the existence. In terms of conservation of energy, all the numerical models are capable of preserving this quantity except [26]. All the numerical models have a consistency of the second order in both space and time, and they are all stable and quadratically convergent in both space and time, for sufficiently small values of the computer parameters. Finally, in terms of implementation, the present numerical methodology is easy to implement computationally, in view that the code only requires solving a couple of linear systems at each time steps. On the other hand, the finite-difference schemes reported in [20,26] are harder to code. This is again due to the nonlinear nature of the numerical models. Since the solution functions are complexvalued, implementations of Newton's method are useless in this point. To overcome this shortcoming, those methods were implemented using a recursive approach.