#### Abstract

In this paper, a new numerical method named Barycentric Lagrange interpolation-based differential quadrature method is implemented to get numerical solution of 1D and 2D coupled nonlinear Schrödinger equations. In the present study, spatial discretization is done with the aid of Barycentric Lagrange interpolation basis function. After that, a reduced system of ordinary differential equations is solved using strong stability, preserving the Runge-Kutta 43 method. In order to check the accuracy of the proposed scheme, we have used the formula of error norm. The matrix stability analysis method is implemented to test the proposed method’s stability, which confirms that the proposed scheme is unconditionally stable. The present scheme produces better results, and it is easy to implement to obtain numerical solutions of a class of partial differential equations.

#### 1. Introduction

##### 1.1. Coupled Schrödinger Equation in One Dimension

Coupled Schrödinger equation in one dimension is presented as follows:

##### 1.2. Coupled Schrödinger Equation in Two Dimensions

Many studies have been done on the notion of nonlinear pulse propagation for more than 30 years. Using the concept of optical solitons in high-speed telecommunication systems was first introduced in 1973 [1, 2]. Later on, when the fiber technique approached some new dimensions, the interest and importance of optical solitons transmission enhanced rapidly. Transmission of multiple gigabytes more than 10,000 km was noticed regularly in 1919. In order to increase the transmission distance, different new approaches like wavelength division multiplexing, distribution management, and polarization division multiplexing came into existence. These newly developed techniques gave rise to the interest of the theoretical modeling related to the propagation of the pulses and the collision of the related systems of communications. In the case of an ideal fiber, the modeling of optical solitons can be represented by the nonlinear Schrödinger system of equations, where these solutions’ behavior is completely known [3, 4], but it is obvious that the real optical fibers have birefringent properties; in such case, the transmission of the pulse happens at a little bit different speed along the two polarization axes. The presence of coupled nonlinear Schrödinger equation also happens in crystals propagation as well as water waves. In [5–7], vector solitons collisions are discussed in coupled nonlinear Schrödinger system of equations and a discussion about collisions with the inelastic property also given by authors, where transmission was represented numerically and explanation of some concepts, is also provided by using basic analytical estimation based on the spatial resonance between two solitons. Yang and Tan [8] observed that the collision represents a fractal structure along with the specific parameters zones; if there will be a little bit of alteration in collision velocity, it can cause a dramatic change in the collision outcome, but in the numerical simulation approach, invariance property was usually neglected in the coupled nonlinear Schrödinger system, like square conservation, playing a vital role in understanding optical fibers’ properties. Since in most of the cases finding exact solution of coupled nonlinear Schrödinger equation is not possible due to its complex nature, most of the researchers have made their contribution to finding a numerical approximation of coupled nonlinear Schrödinger equation. Sun et al. [9] gave the novel six-point approach of finding a numerical approximation of the coupled nonlinear Schrödinger equation. Ismail et al. [10] presented the Alternating Direction Implicit (ADI) method for finding the solution of 2D coupled nonlinear Schrödinger equation. An iterative approach was considered to deal with this equation. Ismail and Taha [11] presented the numerical aspect for approximating coupled nonlinear Schrödinger equation, and, in their paper, finite difference method was implemented to get the solution. Dehghan et al. used [12] meshless local Petrov-Galerkin (MLPG) scheme to approximate the solution of N-coupled nonlinear Schrödinger equation numerically. Xu and Shu [13] presented the local discontinuous Galerkin method for getting the solution of nonlinear Schrödinger and coupled nonlinear Schrödinger equation. In their paper, stability is also fetched for both categories of equations. Aydin and Karasozen [14] presented the approach of multisymplectic integration for the coupled nonlinear Schrödinger equations, where it was mentioned that the six-point scheme is equivalent to the multisymplectic Preissmann approach and the elastic and inelastic properties of solitons were also presented via numerical solutions. Gao et al. [15] gave the implementation of coupled (2 + 1)-dimensional nonlinear Schrödinger equations having variable coefficients related to the optical fibers. Zhang et al. [16] gave the notion of symbolic computation to solve the (2 + 1)-dimensional coupled nonlinear Schrödinger equation and represented the analytical one-soliton and two-soliton solutions under the prescribed conditions using the Hirota method. Wang [17] represented the split-step difference scheme to deal with nonlinear Schrödinger equation as well as coupled nonlinear Schrödinger equation. Ismail and Taha [18] gave a linearly implicit conservative regime to deal with coupled nonlinear Schrödinger equation numerically. The method used by them was accurate in second order in space and time as well. Ismail [19] gave the numerical approximation of coupled nonlinear Schrödinger equation by using the Galerkin approach. In that paper, finite element notion was implemented to deal with the equation. Sonnier and Christov [20] represented the conservative scheme approach for the coupled Schrödinger equation. Ismail [21] gave a method named fourth-order explicit scheme for the coupled nonlinear Schrödinger equation. In that paper, a reduced ODE system of equations was dealt with using the fourth-order Runge-Kutta method, and the boundary conditions used were Neumann and periodic. Sweilam and Al-Bar [22] gave the notion of a variational iteration method for the coupled nonlinear Schrödinger equation. Sun and Qin [23] presented a multisymplectic approach to solving the coupled 1D nonlinear Schrödinger equation. For this purpose, a new six-point scheme almost equivalent to the multisymplectic Preissmann integrator was determined. Abazari and Abazari [24] gave the notion of differential transform method to deal with the numerical aspect of coupled nonlinear Schrödinger equation.

The differential quadrature method (DQM) was first introduced by Bellman et al. [25] in 1972. After Bellman’s work related to DQM, enhancing the depth of this concept goes to the different researchers. Quan and Chang [26, 27] improvised the available approaches of finding weighting coefficients in DQM. The main breakthrough in finding the weighting coefficients required in DQM is the work by Shu and Richards [28]. After that Shu [29] gave the explicit formulation of finding the weighting coefficients. Shu and Chew [30] as well as Shu and Xue [31] gave some more algebraic formulations to determine the weighting coefficients. One more type of DQM is polynomial DQM (PDQM), which is also available in the literature, based upon Lagrange polynomial interpolation. In previous decades, a lot of considerable research has been done by many researchers in this field. Jiwari et al. [32] used weighted average DQM to get the solution of Burgers’ equation numerically when the time derivative was discretized by using the forward difference method and obtained system of linear equations was solved by the Gauss elimination method. Mittal and Dahiya [33] used the notion of MCB-DQM to solve the class of nonlinear viscous wave equations. The MCB spline was used as a basis function in their paper, and the reduced ODE system was solved by RK 43 scheme. The matrix stability analysis method was also implemented in order to verify the stability of the proposed method. Arora and Singh [34] used the concept of MCB-DQM for solving. Burgers’ equation where MCB spline was used in space and SSP-RK43 scheme was implemented to solve a reduced system of ODE. Mittal and Dahiya [35] used the concept of MCB-DQM to solve the hyperbolic diffusion equation numerically, and a reduced set of ODE was dealt with by the SSP-RK43 scheme. Jiwari et al. [36] implemented the notion of polynomial DQM to get the numerical solution of 2D sine-Gordon solitons, where the RK 4 method was implemented to solve the obtained system of ODE. Jiwari et al. [37] used the DQM to solve the 1D hyperbolic telegraph equation of second order numerically. In their paper, polynomial DQM was implemented, and obtained system of ODE was treated with the RK 4 method. Jiwari et al. [38] used the notion of polynomial DQM to get the numerical solution of linear hyperbolic telegraph equation of 2 dimensions with Dirichlet and Neumann boundary conditions. Obtained system of ODE was solved by using the RK 4 method. Mittal and Bhatia [39] used a modified B-spline in DQM to solve the 2D hyperbolic telegraph equation. In their paper, the resulting system of ODE was solved by the SSP-RK43 method. Shukla et al. [40] used the notion of exponential modified cubic B-spline basis function in DQM to get numerical solution of 3D nonlinear wave equation. The obtained system of ODE was solved by using the SSP-RK54 scheme. The stability of the method was investigated by the matrix stability analysis method. Korkmaz and Dag [41] used the notion of DQM to solve the nonlinear Schrödinger equation. In their paper, the differential quadrature method was developed by using cosine expansion. Korkmaz and Dag [42] implemented the notion of DQM to get numerical solutions of complex modified Korteweg-de Vries equation. Korkmaz and Dag [43] implemented the notion of cubic B-spline-based DQM to get a numerical solution of the advection-diffusion equation. Korkmaz and Dag [44] implemented the concept of Sinc DQM to solve shock wave simulations. In their research paper, Sinc DQM was used to discretize spatial variables, and the four-stage RK method was used for time discretization. Mittal and Rohila [45] used modified cubic DQM to get a numerical solution of the reaction-diffusion system. Mittal and Dahiya [46] used the concept of cubic B-spline in DQM to get the solution of the 3D telegraphic equation.

The Barycentric Lagrange interpolation formula is obtained from the well-known Lagrange interpolation formula. Barycentric form of Lagrange polynomial interpolation was proposed by Berrut et al. [47]. Berrut also gave the review of Barycentric rational interpolation and its numerous applications in [48]. Guttel presented a convergent scheme of Barycentric rational interpolation formula in [49]. In various cases, the Lagrange interpolation formula is a good choice to tackle polynomial interpolation. The Lagrange interpolation formula can be modified with the aid of Barycentric interpolation. Lagrange interpolation is one of the most commonly used approaches for the analysis of theory purpose. Still, it is usually avoided from the numerical perspective as it is complex in computing and numerical instability. A Barycentric form of Lagrange polynomial interpolation can reduce the complexities presented in the original Lagrange interpolation formula and give rise to interpolation accuracy. The Barycentric form of Lagrange polynomial interpolation function is related to the notion of the continuous function , *x* [−1, 1], at the given (*n* + 1) distinct nodes , *i* = 0, 1, …, *n*, and it is constructed as follows:wherewhere is denoted as the modified formula for the Lagrange polynomial interpolation function. are the functional values at the given interpolation nodes . are Barycentric Lagrange basis functions, which follow the following property:where is known as the Kronecker-delta function and ; and are known as the corresponding Barycentric weights:

A lot of work has been reported in the literature regarding the Barycentric form of the Lagrange interpolation formula, and numerous techniques are presented for the numerical solution of partial differential equations. Ma et al. [50] reported the concepts of a meshless collocation method based upon Barycentric rational interpolation for hyperbolic telegraph equation in 2D. In this paper, Barycentric rational interpolation was used for spatial discretization, and a system of second-order ODE was obtained based upon the collocation method. Liu et al. [51] reported that the collective compact theory showed the notion of the Barycentric interpolation-based collocation method for the solution of linear and nonlinear high-dimensional Fredholm integral equation as well as the convergence of the method. Liu et al. [52] implemented a mesh-free aspect by using the Barycentric Lagrange interpolation formula to solve the multidimensional system of Fredholm integral equations. This used method was an improvised form of Lagrange polynomial interpolation, which has higher accuracy. Berrut and Klein [48] presented recent development regarding linear Barycentric rational interpolation. Higham [53] presented numerical stability for the Barycentric form of Lagrange polynomial interpolation. Berrut and Mittelman [54] gave the notion of different matrices for directly determining Barycentric weights regarding rational interpolation. Berrut and Trefethen [47] discussed Barycentric Lagrange interpolation and stated that it is much faster than the original one. Lawrence and Corless [55] discussed the notion of stability of root-finding regarding Barycentric Lagrange polynomial interpolation. Luo et al. [56] implemented the Barycentric rational collocation scheme to solve a class of parabolic PDEs, nonlinear in nature. Floater and Hormann [57] discussed Barycentric form-based rational interpolation having no poles and higher rates on approximation. Wu et al. [58] implemented the notion of a Barycentric polynomial interpolation-based collocation scheme to get the numerical solution of a class of nonlinear PDEs.

In the present paper, the Barycentric form of Lagrange interpolation-based differential quadrature method is applied to solve the 1D and 2D coupled nonlinear Schrödinger equations. By reviewing the literature in depth, the authors found that this notion has never been implemented to solve the 1D and 2D coupled nonlinear Schrödinger equations. The present paper is divided into different sections. In Section 2, a numerical scheme is presented to approximate the solution of 1D and 2D coupled nonlinear Schrödinger equations. A generalized formula is developed to find the weighting coefficients of the differential quadrature method. We have considered the nonuniform mesh in the present paper by taking Chebyshev–Gauss–Lobatto grid points, as, at these grid points, better results are obtained. In Section 3, four numerical examples are discussed. In Section 4, the present scheme’s stability is discussed using the matrix stability analysis method.

#### 2. Numerical Scheme (Barycentric Lagrange Polynomial Interpolation DQM)

The Barycentric formula for a continuous function for the domain [−1, 1] at the node points in the interval [−1, 1] is given as follows [51]:where is known as the modified form of Lagrange polynomial interpolation, are corresponding functional values at given interpolation node points, and mentioned factors are known as Barycentric Lagrange basis functions:and property is also satisfied, where is the associated Kronecker-delta function and . In order to develop a generalized formula, the method is explained with the number of node points. In every case at different node points, different basis functions are applied in the formula of DQM to obtain the corresponding weighting coefficients so that the partial derivative can be approximated.

In the case of grid points, , there will exist basis functions, , and weighting coefficients will exist.

Derivative approximation of function in DQM is given by the following formula:

Using different Barycentric basis functions along with provided weights in the formula of DQM, the following generalized formula will be achieved. Different weighting coefficients used in the differential quadrature method can be fetched.

In the present paper, DQM is implemented to approximate the spatial derivative using the Barycentric Lagrange interpolation basis function.where are the nondiagonal weighting coefficients when and are the coefficients when . In first-order derivative approximation can be obtained from the generalized formula of Barycentric Lagrange interpolation.

In the approximation of -order partial derivative, used weighting coefficients can be obtained by using the following formula [59]:where is weighting coefficient of the first-order partial derivative and is weighting coefficient of the order derivative. Weighting coefficient in second-order partial derivative approximation can be obtained by using the following formula [59]:where is the weighting coefficient of second-order partial derivative approximation.

In this complete numerical study of the present paper, we have given preference to the Chebyshev–Gauss–Lobatto grid points in order to find the numerical solution of coupled 1D and coupled 2D Burgers’ equation. The mentioned Chebyshev–Gauss–Lobatto grid points are given as follows:

#### 3. Numerical Examples and Discussion

In the present section, we have discussed four numerical experiments. Among these four numerical experiments, the first numerical example is concerned with 1D coupled nonlinear Schrödinger equation and the other three examples are related to 2D coupled nonlinear Schrödinger equation. In Table 1, error norms are represented for at different time levels and for Δ*t* = 0.01, 0.05, and 0.08, respectively. In Table 2, error norm is calculated at different time levels and for Δ*t* = 0.01, 0.05, and 0.08, respectively. In Table 3, a comparison is shown between error norms with [11]. The present results are much better than the compared ones. In Figure 1, numerical plot is given for first wave amplitude in the form of surface and contour plots for *N* = 41, , , , , and at time levels *t* = 1 and 2, respectively. In Figure 2, numerical plot is given for first wave amplitude in the form of surface and contour plots for *N* = 41, , , , , and at time levels *t* = 3 and 4, respectively. In Figure 3, numerical plot is given for second wave amplitude in the form of surface and contour plots for *N* = 41, , , , , and at time levels *t* = 1 and 2, respectively. In Figure 4, numerical plot is given for second wave amplitude in the form of surface and contour plots for *N* = 41, , , , , and at time levels *t* = 3 and 4, respectively. In Figure 5, surface and contour plots are provided for first wave amplitude for *N* = 11, , and at time levels 1 and 2, respectively. In Figure 6, surface and contour plots are provided for first wave amplitude for *N* = 11, , and at time levels 3 and 4, respectively. In Figure 7, surface and contour plots are provided for second wave amplitude for *N* = 11, , and at time levels 1 and 2, respectively. In Figure 8, surface and contour plots are provided for second wave amplitude for *N* = 11, , and at time levels 3 and 4, respectively. In Figure 9, surface and contour plots are provided for first wave amplitude for *N* = 11, , and at time levels 1 and 2, respectively. In Figure 10, surface and contour plots are provided for first wave amplitude for *N* = 11, , and at time levels 3 and 4, respectively. In Figure 11, surface and contour plots are provided for second wave amplitude for *N* = 11, , and at time levels 1 and 2, respectively. In Figure 12, surface and contour plots are provided for second wave amplitude for *N* = 11, , and at time levels 3 and 4, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

*Example 1. *Single soliton:

In this example, 1D coupled nonlinear Schrödinger equations (1) and (2) are considered with the following initial conditions:where *e* and both are considered as constant values. Boundary conditions imposed are natural boundary conditions:

*Example 2. *In this example, 2D coupled nonlinear Schrödinger equations (3) and (4) are given with following initial conditions:where computational domain = [−10, 10], .

*Example 3. *In this example, 2D coupled nonlinear Schrödinger equations (3) and (4) are considered with the following initial conditions [62]:where boundary conditions arewhere consideration is = = −1 and = = −1. Computational domain of this problem is [−1, 1] [−1, 1].

*Example 4. *In this example, 2D coupled nonlinear Schrödinger equations (3) and (4) are considered with the following initial conditions [62]:where boundary conditions arewhere consideration is = = −1 and = = −1. Domain of computation of this problem is [−1, 1] [−1, 1]. In this example, no exact solution is provided. The numerical solution, in this case, is tried to depict with the aid of graphs.

#### 4. Stability

The following system of equations will be obtained by implementing the discretization formulae of the partial derivatives in equations (1) and (2):

The above system of equations can be written as follows:where a matrix is fetched from the abovementioned system of ordinary differential equations, and is the corresponding nonlinear term. The stability of the above system depends upon the eigenvalues of matrix .where

By implementing the discretization formula in equations (3) and (4), a newly formed system of ordinary differential equations will be obtained as follows:where the above system of ordinary differential equations can be written as follows:

The stability of the present scheme will depend upon the eigenvalues of matrix and is the corresponding nonlinear term.where

The concept of the proposed method’s stability is discussed here by taking the basic concepts from [60, 61]. The obtained systems of ordinary differential equations will be stable if the eigenvalues of matrices and will lie in the following range:(a)Real : −2.78 < < 0(b)Pure imaginary : −2 < < (c)Complex : will lie inside the region as per literature [61].

The stability of the proposed scheme is shown with the aid of Figures 13 and 14 at the range of the number of grid points. As per the stability criteria, it is quite clear that the developed scheme is unconditionally stable.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

#### 5. Conclusion

In the present paper, the Barycentric Lagrange interpolation-based differential quadrature method is applied to get numerical solution of 1D and 2D coupled nonlinear Schrödinger equations. In most cases, obtaining analytical solutions of such partial differential equations is a cumbersome task, so finding the numerical solution by novel techniques arises. In the present paper, we have discussed the proposed scheme’s accuracy by taking four numerical experiments. By using the matrix stability analysis method, it is confirmed that the present scheme is unconditionally stable. By all these observations, it can be concluded that the proposed scheme produces far better results than existing results. This scheme will help the researchers to explore some new dimensions to find numerical solutions to a class of partial differential equations.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest to report regarding the present study.

#### Acknowledgments

The authors would like to acknowledge the support of Taif University Researchers Supporting Project (no. TURSP-2020/211), Taif University, Taif, Saudi Arabia.