# Optimized integrating factor technique for Schrödinger-like equations

The integrating factor technique is widely used to solve numerically (in particular) the Schrödinger equation in the context of spectral methods. Here, we present an improvement of this method exploiting the freedom provided by the gauge condition of the potential. Optimal gauge conditions are derived considering the equation and the temporal numerical resolution with an adaptive embedded scheme of arbitrary order. We illustrate this approach with the nonlinear Schrödinger (NLS) and with the Schrödinger-Newton (SN) equations. We show that this optimization increases significantly the overall computational speed, sometimes by a factor five or more. This gain is crucial for long time simulations.

## Authors

• 1 publication
• 1 publication
• 1 publication
10/23/2020

### Maximum bound principle preserving integrating factor Runge-Kutta methods for semilinear parabolic equations

A large class of semilinear parabolic equations satisfy the maximum boun...
09/09/2019

### A comparison between pre-Newton and post-Newton approaches for solving a physical singular second-order boundary problem in the semi-infinite interval

In this paper, two numerical approaches based on the Newton iteration me...
04/20/2021

### A decoupling and linearizing discretization for poroelasticity with nonlinear permeability

We analyze a semi-explicit time discretization scheme of first order for...
04/20/2021

### Stabilized integrating factor Runge-Kutta method and unconditional preservation of maximum bound principle

Maximum bound principle (MBP) is an important property for a large class...
04/07/2020

### Policy iteration for Hamilton-Jacobi-Bellman equations with control constraints

Policy iteration is a widely used technique to solve the Hamilton Jacobi...
09/30/2021

### A Fast Robust Numerical Continuation Solver to a Two-Dimensional Spectral Estimation Problem

This paper presents a fast algorithm to solve a spectral estimation prob...
09/30/2020

### A frequency-dependent p-adaptive technique for spectral methods

When using spectral methods, a question arises as how to determine the e...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

The Schrödinger equation is used in many field of Physics. In dimensionless units, it takes the form

 i∂tψ + 12∇2ψ − Vψ = 0, (1)

where is a function of the spatial coordinates and of the time , is the Laplace operator, and the potential is generally a function of space and time and, possibly, a functional of . In this paper, we focus more specifically on the nonlinear Schrödinger equation ( proportional to ) and on the Schrödinger–Newton (or Schrödinger–Poisson) equation in which the Laplacian of the potential is proportional to . These special cases were chosen for clarity and because they are of practical interest, but the method presented here can be extended to more general potentials (and equations).

With the nonlinear Schrödinger equation (NLS) considered here, the (local nonlinear) potential is

 V = g|ψ|2, (2)

where is a coupling constant. For the interaction is repulsive, while, for it is attractive. The NLS describes various physical phenomena, such as Bose–Einstein condensates [1], laser beams in some nonlinear media [2], water wave packets [3], etc.

With the Schrödinger–Poisson equation, the potential is given by the Poisson equation

 ∇2V = g|ψ|2, (3)

where is another coupling constant, the interaction being attractive if and repulsive if . It is therefore nonlinear and non-local, giving rise to collective phenomena [4], appearing for instance in optics [5, 6, 7], Bose–Einstein condensates [8], cosmology [9, 10, 11] and theories describing the quantum collapse of the wave function [12, 13]. It is also used as a model to perform cosmological simulations in the semi-classical limit [14].

The above equations cannot be solved analytically (except for very special cases) and numerical methods must be employed. In this paper, we focus on spectral methods for the spatial resolution, i.e., methods that are based on fast Fourier transform (FFT) techniques, that are specially efficient and accurate

[15]. For the temporal resolution, two families of methods are commonly employed to solve Schrödinger-like equations: integrating factors [16] and split-step integrators [17]. The latter methods have been used to integrate both the SN and NLS equations, but the former is used essentially to solve the NLS, with very performing results [18, 19]. In this note, we focus on the former technique, which consists in integrating analytically the linear part of the equation and integrating numerically the remaining nonlinear part with a classical method [16]. The principle of the method is described as follow.

Writing the Schrödinger equation in the generic form

 i∂tψ = F(r,t,ψ),ψ = ψ(r,t), (4)

the right-hand side is split into linear and nonlinear parts

 i∂tψ + Lψ = N(r,t,ψ), (5)

where is an easily computable autonomous linear operator and is the remaining (usually) nonlinear part. At the -th time-step, with , considering the change of dependent variable

 ϕ \tiny{def}= exp[(t−tn)L]ψ⟹i∂tψ = exp[(tn−t)L](i∂tϕ−Lϕ), (6)

so at , the equation (5) is rewritten

 i∂tϕ =exp[(t−tn)L]N. (7)

The operator being well chosen, the stiffness of (5) is considerably reduced and the equation (7) is (hopefully) well approximated by algebraic polynomials for . Thus, standard time-stepping methods, such as an adaptive Runge–Kutta method [20, 21], can be used to efficiently solve (7

). To do so, the solution is evaluated at two different orders and a local error is estimated as the difference between those quantities. Popular integrators can be found in

[20, 22].

It is straightforward to apply this strategy to the Schrödinger equation (1) since and the potential are, respectively, linear and nonlinear operators of . By switching to Fourier space in position, the equation becomes

 i∂tˆψ − 12k2ˆψ − ˆ\/Vψ\/ = 0, (8)

where hats denote the Fourier transform of the underneath quantity and (

the wave vector). The equation is now in a form where the application of the integrating factor (IF) method is straightforward, i.e., (

8) becomes

 i∂tϕ = −\/iexp[i2\/k2\/(t−tn)\/]ˆ\/Vψ\/, (9)

where . If the nonlinear part of the equation is zero, then and any (reasonable) temporal scheme will produce the exact solution . In other words, the integrating factor technique is exact for linear equations. This indicates that the numerical errors depend on the magnitude of the nonlinear part. Therefore, in order to minimise these errors, a strategy consists in minimizing the magnitude of at each time-step. To do so, we exploit the gauge invariance of the Schrödinger equation: if is a solution of (1) at a given time , then is a solution of

 i∂tΨ + 12∇2Ψ − (V+C)Ψ = 0, (10)

as one can easily verify. Thus, at each time-step, adding a constant to in (1), modifies the solution as

 ψ(tn) → ψ(tn)e−\/i\/φ,φ \tiny{def}= n∑j=1Cjhj, (11)

where is the -th time-step. Of course, at the end of the computations, the operation (11) can be easily reverted if the original phase is relevant. Using this procedure, we observed up to a five-fold speed increase (the overall computational time is divided by about five) compared to taking . Of course, the speed-up varies depending on the initial condition, of the (spatial and temporal) numerical schemes and on the choice of gauge corresponding to .

In this paper, we derive some analytic formulas giving an optimal in order to maximise the time-step, i.e., to minimize the overall computational time of the numerical resolution. We emphasize that the thus obtained optimal values of do not affect the accuracy of the numerical solution, leaving it unchanged with respect to the case. Two strategies are presented. In section 2, a first ‘natural’ approach to derive a suitable is based on the analytical structure of the equation and it is independent of the numerical algorithm employed for its resolution. More precisely, is obtained minimizing a norm of the right-hand side of the equation (5). This provides an easy manner to obtain a formula that is moreover computationally ‘cheap’. This expression is however only near optimal, so a ‘better’ expression is subsequently derived. Considering both the equations and the numerical algorithms, a second optimal expression for is derived in section 3. This approach consists in minimizing exactly the numerical error and thus explicitly dependents on the numerical scheme. This provides a more accurate, but computationally expensive, solution. The advances of these special choices are illustrated numerically in section 4. Finally, a summary and perspectives are drawn in section 5.

## 2. Near optimal Cn

As mentioned above, if properly chosen, the integrating factor is able to reduce the stiffness of the equation, making the numerical integration more efficient. In addition, the magnitude of the nonlinear part of (7) also contributes to the efficiency of the numerical integration. Specifically, if is zero, and the integrating factor technique is exact. Thus, the efficiency of the algorithm is expected to increase as the magnitude of gets smaller, and subsequently the overall computational time should be reduced. Here, we show how to choose the arbitrary constant in order to reduce the magnitude of the nonlinear part . In the case of the Schrödinger equation, we have

 N(k,t;ϕ;Cn) = −\/iexp[i2\/k2\/(t−tn)\/]F{(V+Cn)ψ}, (12)

where denotes the Fourier transform and .

A natural strategy is to minimize the -norm, namely

 Gn(Cn)% \tiny{def}=1M[M/2]−1∑m=−[M/2]|N(km,tn;ϕ;Cn)|2, (13)

where is the number of spatial modes, square brackets denote the integer part and is the -th Fourier mode. The explicit expression of can be found exploiting the definition of the discrete Fourier transform. For simplicity, we do the calculations in one dimension (1D) without loss of generality, since the final result is independent of the spatial dimension . From Parseval theorem, one obtains

 Gn(Cn) =[M/2]−1∑ℓ=−[M/2](Vℓ+Cn)2|ψℓ|2, (14)

where and at time . Since the function is a second-order polynomial in , it admits an unique minimum, which is obtained from the equation , yielding

 Cn=−⎛⎝[M/2]−1∑ℓ=−[M/2]Vℓ|ψℓ|2⎞⎠/⎛⎝[M/2]−1∑ℓ=−[M/2]|ψℓ|2⎞⎠\tiny{def}=~Cn. (15)

Therefore, at each time step , , which is the value of minimizing the -norm of , is obtained from (15). We show below that even though this approach is not unique (i.e., different norms could be considered), the provided solution is quite advantageous compared to others, being computationally cheap and independent on the order of the numerical scheme.

## 3. Optimal Cn

We show here another way to choose the arbitrary constant in order to improve the algorithm efficiency and reduce the overall computational time. This approach is based on the principles of the adaptive time-step procedure, where at each time step , an error between two approximated solutions of different orders is estimated. Since the smaller this quantity the larger the time-step, minimizing allows to choose a larger time-step, speeding-up the numerical integration and keeping roughly the same numerical error. More specifically, the error depends on the arbitrary constant , hence the minimization can be performed (see below) choosing an appropriate . In this section, we first recall the method for determining the size of the time step used in the Runge–Kutta procedures; interested readers should refer to [20] for further details. Although the determination of can be formally presented for any embedded Runge–Kutta schemes, this results in very cumbersome calculations with little insights. Thus, for brevity and clarity, we illustrate the method with the Heun method (that is a second-order Runge–Kutta method with an embedded first-order explicit Euler scheme for the time stepping [23]). We then sketch-out how this procedure can be implemented for generic embedded Runge–Kutta methods.

### 3.1. Principle of the adaptive time-step procedure

For the time stepping, embedded Runge–Kutta methods estimate the quadrature error comparing the results of two orders of the time integrator [20]. For a solver of order with an embedded -order scheme (hereafter schemes of orders ), at the -th time step, the error is [24]

 Δn \tiny{def}=    ⎷1M[M/2]−1∑m=−[M/2]⎛⎜ ⎜⎝∣∣ϕ(km,tn)−˜ϕ(km,tn)∣∣\sc Tol+max(|ϕ(km,tn)|,∣∣~ϕ(km,tn)∣∣)×% \sc Tol⎞⎟ ⎟⎠2, (16)

where is the number of spatial modes, square brackets denote the integer part, is the -th order solution at the -th Fourier mode, the tilde notation indicating the solution at order , and Tol is the tolerance (parameter defining the desired precision of the time-integration). The time step is accepted if the error is smaller than the tolerance Tol, otherwise is reduced and this step is recomputed. being accepted, the next time step is obtained assuming the largest error equal to the tolerance. In order to avoid an excess of rejected time steps, we use the Proportional Integral (PI) Step Control [24], which chooses the optimal time step as

 hn+1 = hnΔ−bnΔcn−1, (17)

where , , being the order of the chosen integrator [25]. Interested readers should refer to [20] for details on this classical procedure.

### 3.2. Optimum time step

Since the constant can be chosen freely, we seek for the value of providing the largest , namely, to maximize the right-hand side of (17). Since and are determined at the previous time-step, only in (17) depends on . Thus, in order to maximize , must be minimized, i.e., one must solve . This derivation being characterized by cumbersome algebra for general embedded Runge–Kutta schemes, we illustrate the case of the Heun algorithm (that is a second-order Runge–Kutta method with an embedded first-order explicit Euler scheme for the time stepping [23]), the principle being the same for higher order integrators. Also for simplicity, we give the calculations in one dimension (1D) without loss of generality, since the final result is independent of the spatial dimension.

#### 3.2.1. Optimum Cn for Heun’s method

Heun’s method consists, here, in solving the initial value problem (for )

 i∂tϕ = f(k,t;ϕ;Cn)\tiny{def}=−\/iexp[i2k2(t−tn)]F{(V+Cn)ψ}, (18)

and

 ϕ(k,t) \tiny{def}=exp[i2k2(t−tn)]F{ψ(x,t)}. (19)

Hereafter, for brevity, we denote

 ϕn=ϕn(k) \tiny{def}= ϕ(k,tn),ψn=ψn(x)\tiny{def}= ψ(x,tn),Vn=Vn(x)\tiny{def}=V(x,tn). (20)

At time , the first- and second-order (in ) approximations of , respectively and , are

 ˜ϕn+1 (21) ϕn+1 = ϕn + 12hn[f(k,tn;ϕn;Cn)+f(k,tn+hn;ϕn+hn\/f(k,tn;ϕn;Cn);Cn)]. (22)

The next time-step is chosen using equation (17). For our equation, the difference between the first- and second-order approximations is such that

 (Δϕn+1)2 = 14h2n∣∣f(k,tn;ϕn;Cn)−f(k,tn+hn;ϕn+hn\/f(k,tn;ϕn;Cn);Cn)∣∣2 = 14h2n∣∣f(k,tn;ϕn;Cn)+iei\/k2\/hn/2× F{(Vn+1+Cn)\/F−1{e−i\/k2\/hn/2(ϕn+hnf(k,tn;ϕn;Cn))}}\/∣∣2, (23)

where . We note that the absolute value in (3.2.1) is of first-order in , as one can easily check with a Taylor expansion around , so . More precisely, after some elementary algebra, one finds

 (Δϕn+1)2 = 14h4n∣∣F{(Vn+Cn)2ψn} + iF{∂tVnψn} + F{∂xVn∂xψn} (24) + 12F{∂xxVnψn} ∣∣2 + O(h5n),

which, defining

 α(x,t;Cn) \tiny{def}=(Vn+Cn)2ψn,β(x,t)\tiny{def}=i∂tVnψn + ∂xVn∂xψn+ 12∂xxVnψn, (25)

can be rewritten as

 (Δϕn+1)2 = 14h4n|F{α(Cn)+β}|2 + O(h5n). (26)

 En(Cn)\tiny{def}= 1M[M/2]−1∑m=−[M/2]Δϕ2n+1(km,tn;Cn), (27)

substituting (26) into (27) and exploiting the definition of the discrete Fourier transform, one obtains (using Parseval theorem)

 En(Cn) =1M[M/2]−1∑m=−[M/2]∣∣ ∣∣[M/2]−1∑ℓ=−[M/2]e−\/2\/i\/π\/m\/ℓ/M(αℓ(Cn)+βℓ)∣∣ ∣∣214h4n + O(h5n) =[M/2]−1∑ℓ=−[M/2]|αℓ(Cn)+βℓ|214h4n+ O(h5n). (28)

The minimum of , obtained from the equation , is such that

 [M/2]−1∑ℓ=−[M/2]d|αℓ(Cn)|2dCn + 2Re(dαℓ(Cn)dCnβ∗ℓ)=0, (29)

Therefore, the optimum providing the largest , in the case of Heun’s method, is a solution of (29).

#### 3.2.2. Optimum Cn for generic embedded Runge–Kutta schemes

The optimum for general embedded Runge–Kutta schemes can be obtained following the same principles illustrated above with the Heun algorithm. However, the algebraic calculations get rapidly very cumbersome, leading to expensive computations that, in most cases, exceeds the time gained with a larger step. Here, we sketch-out the procedure for generic embedded Runge–Kutta methods, considering solvers of order with an embedded -order scheme (for other embedded or extrapolation methods, the procedure is completely analogue). For a -stage method, the error can be written as [24]

 (Δϕn+1)2=∣∣ ∣∣s∑ℓ=1dℓwℓ∣∣ ∣∣2,dℓ\tiny{% def}=as,ℓ−bℓ,wℓ% \tiny{def}=hnf(k,tn+cℓhn;ϕn+ℓ−1∑r=1aℓ,rwr;Cn), (30)

where , and are the coefficients of the Butcher tableau which characterizes the integrator [24]. Using Taylor expansions and un-nesting the scheme, it is possible to prove that a result with a similar structure compared with (26) is obtained. In this case, the number of stages appears as exponent in the function , which takes the form . The function , on the other hand, becomes explicitly dependent on , involving a number of terms growing exponentially with . For this reason, even though the exact result can always be achieved, the computational time needed to minimize the error (16) is often larger than the time gained with a larger step, especially for higher order schemes (). In the next section, we show how, for practical applications, an exact solution is not necessary to improve the algorithm, and (15) represents a fast and accurate method.

## 4. Numerical examples

Here, we consider numerical examples where we apply this method, focusing on both the SN and NLS equations solved with the Dormand and Prince 5(4) integrator [22] in one and two spatial dimensions. In all cases, we set open boundary conditions for the potential, while the initial conditions and the value of the physical parameters are chosen to be very close to regimes of physical interest, as described in [1, 2, 3]. Details on the numerical simulations can be found in A. The gain factor provided by the method depends on the optimal value of compared to the case, which changes from case to case as a function of the boundary conditions for the potential and of the profile of the solution. Specifically, since the gain factor is evaluated with respect to the case, the more the optimal value of is far from zero, the larger the gain factor gets.

For the one-dimensional NLS, some analytical stationary solutions are known. We then use one of these solutions (see A) as initial condition. For all other cases (SN and NLS 2D), no such stationary solutions are known, so we use gaussian initial conditions.

In Figure 1, we show the average time-step , for an entire simulation with time steps, as a function of for the one-dimensional SN and NLS equations. These plots are generated taking constant for the entire simulations, in order to better appreciate the strong dependence of the time-step on the choice of the gauge for the potential. In Figure 2, we report the result of simulations performed choosing the near optimal at each time-step. Note that, for the one-dimensional NLS, the solution being stationary, and the optimum do not change in time, that is not the case in 2D. We show the time-step as a function of time for the one-dimensional SN and NLS equations, comparing the case with . In both cases, the time-step chosen by the algorithm with the optimization of the gauge constant proves to be larger, compared to the case. In Table 1, we show the number of time-loops required to run each simulation and the time needed to run the simulation (in seconds) for the cases and . For NLS, in the one dimensional case we achieve roughly a improvement in terms of speed gain between the and cases, while in two dimensions the speed gain is only approximately since, here, the value of is very close to zero. As long as the equation is concerned, (15) proved to reduce remarkably both the number of time loops and the effective time for the total simulation, providing up to a factor of improvement with respect to the case in 1D and up to a factor in 2D.

## 5. Conclusion

Exploiting a gauge condition on the potential, we optimized the integrating factor technique applied to the nonlinear Schrödinger and Schrödinger–Newton equations. Although the exact values of the piece-wise constant minimizing the error (16) (therefore maximizing the time-step) is in principle always possible to compute (e.g., with a computer algebra system), its expression depends on the particular numerical scheme chosen and it becomes complicated as the order of the method increases, resulting in a high computational cost. However, the near-optimal value obtained from the first approach we described, based on the minimization of the -norm of the nonlinear part of the equation, proved to be an accurate and efficient solution in the tested cases. Thus, being computationally extremely cheap and independent of the particular numerical scheme employed, this is the approach one should choose for most simulations, at least when the computation of is not very expensive. For Schrödinger-like equations with hard to compute potentials, most of the computational time is spent in the calculation of . For these very demanding equations, the extra cost needed to compute the optimum (instead of the near optimum ) is negligible in comparison, so could be preferable.

For the cases tested here, we found a speed-up in the computation time up to a factor , the speed-up depending on the equation and on the physical regime. These examples show that this approach provides significant speed improvements, that with minor modifications of the original algorithm. Though we focused on the nonlinear Schrödinger and Schrödinger–Newton equations, the method principle is independent on the particular potential considered, so this approach can be extended to other Schrödinger-like equations. More generally, the idea behind the method presented in this note can be, at least in principle, generalized and extended to other equations with similar gauge conditions.

## Appendix A Numerical simulations

For the one-dimensional NLS, we considered the case (see 3) and we used as initial condition. We discretised space with points, in a computational box of length . The two-dimensional NLS, which is often employed in optics to model self-focusing beams in a medium with a cubic non-linearity [26, 27], presents a finite time (blow-up) singularity [28]. More specifically, whenever the initial condition satisfies , the norm of the solution, or of one of its derivatives, becomes unbounded in finite time. For this reason, we stop the simulation at , i.e., before the singularity occurs. We set as initial condition and we consider the case, for which the corresponding initial energy is , hence quite close to the singular regime; for the spatial discretization we used and (squared box with side discretised with nodes). For both the one and two dimensional SN equations, we set and considered a Gaussian initial condition, where is the normalization factor, fixed such that . The parameters of the spatial discretization are and in 1D, while for the 2D case we set and . To solve the SN and the NLS equation we used the Dormand and Prince 5(4) integrator [22].

## References

• [1] Franco Dalfovo, Stefano Giorgini, Lev P Pitaevskii, and Sandro Stringari. Theory of Bose–Einstein condensation in trapped gases. Reviews of Modern Physics, 71(3):463, 1999.
• [2] Akira Hasegawa. Optical solitons in fibers. In Optical Solitons in Fibers, pages 1–74. Springer, 1989.
• [3] Christian Kharif, Efim Pelinovsky, and Alexey Slunyaev. Rogue waves in the ocean. Springer Science & Business Media, 2008.
• [4] James Binney and Scott Tremaine. Galactic Dynamics: Second Edition. 2008.
• [5] FW Dabby and JR Whinnery. Thermal self-focusing of laser beams in lead glasses. Applied Physics Letters, 13(8):284–286, 1968.
• [6] Carmel Rotschild, Oren Cohen, Ofer Manela, Mordechai Segev, and Tal Carmon. Solitons in nonlinear media with an infinite range of nonlocality: first observation of coherent elliptic solitons and of vortex-ring solitons. Physical review letters, 95(21):213904, 2005.
• [7] Carmel Rotschild, Barak Alfassi, Oren Cohen, and Mordechai Segev. Long-range interactions between optical solitons. Nature Physics, 2(11):769–774, 2006.
• [8] F Siddhartha Guzman and L Arturo Urena-Lopez. Gravitational cooling of self-gravitating Bose condensates. The Astrophysical Journal, 645(2):814, 2006.
• [9] Wayne Hu, Rennan Barkana, and Andrei Gruzinov. Fuzzy cold dark matter: the wave properties of ultralight particles. Physical Review Letters, 85(6):1158, 2000.
• [10] Angel Paredes and Humberto Michinel. Interference of dark matter solitons and galactic offsets. Physics of the Dark Universe, 12:50–55, 2016.
• [11] David JE Marsh and Jens C Niemeyer. Strong constraints on fuzzy dark matter from ultrafaint dwarf galaxy Eridanus ii. Physical review letters, 123(5):051103, 2019.
• [12] Lajos Diosi. Gravitation and quantumm-echanical localization of macroobjects. arXiv preprint arXiv:1412.0201, 2014.
• [13] Roger Penrose. On gravity’s role in quantum state reduction. General relativity and gravitation, 28(5):581–600, 1996.
• [14] Lawrence M Widrow and Nick Kaiser. Using the Schrödinger equation to simulate collisionless matter. The Astrophysical Journal, 416:L71, 1993.
• [15] C. Canuto, M. Y. Hussaini, A. Quarteroni, and Th. A. Zang. Spectral Methods. Fundamentals in Single Domains. Scientific Computation. Springer, 2nd edition, 2008.
• [16] J Douglas Lawson. Generalized Runge–Kutta processes for stable systems with large Lipschitz constants. SIAM J. Numer. Anal., 4(3):372–380, 1967.
• [17] Sergio Blanes, Fernando Casas, and Ander Murua. Splitting and composition methods in the numerical integration of differential equations. arXiv preprint arXiv:0812.0377, 2008.
• [18] Philipp Bader, Sergio Blanes, Fernando Casas, and Mechthild Thalhammer. Efficient time integration methods for Gross–Pitaevskii equations with rotation term. arXiv preprint arXiv:1910.12097, 2019.
• [19] Christophe Besse, Geneviève Dujardin, and Ingrid Lacroix-Violet. High order exponential integrators for nonlinear Schrödinger equations with application to rotating Bose–Einstein condensates. SIAM J. Numer. Anal., 55(3):1387–1411, 2017.
• [20] Roger Alexander.

Solving ordinary differential equations i: Nonstiff problems (e. hairer, sp norsett, and g. wanner).

Siam Review, 32(3):485, 1990.
• [21] John Charles Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, 2016.
• [22] John R Dormand and Peter J Prince. A family of embedded Runge–Kutta formulae. J. Comp. Appl. Math., 6(1):19–26, 1980.
• [23] Endre Süli and David F Mayers. An introduction to numerical analysis. Cambridge university press, 2003.
• [24] Gerhard Wanner and Ernst Hairer. Solving ordinary differential equations II. Springer Berlin Heidelberg, 1996.
• [25] Kjell Gustafsson. Control theoretic techniques for stepsize selection in explicit Runge–Kutta methods. ACM Transactions on Mathematical Software (TOMS), 17(4):533–554, 1991.
• [26] VE Zakharov and VS Synakh. The nature of the self-focusing singularity. Zh. Eksp. Teor. Fiz, 68:940–947, 1975.
• [27] Kimiaki Konno and Hiromitsu Suzuki. Self-focussing of laser beam in nonlinear media. Physica Scripta, 20(3-4):382, 1979.
• [28] Catherine Sulem and Pierre-Louis Sulem. The nonlinear Schrödinger equation: self-focusing and wave collapse, volume 139. Springer Science & Business Media, 2007.