A Study of Traveling Wave Structures and Numerical Investigations into the Coupled Nonlinear Schrödinger Equation Using Advanced Mathematical Techniques

Abstract

This article explores adapted mathematical methods to solve the coupled nonlinear Schrödinger (C-NLS) equation through analytical and numerical methods. To obtain exact solutions for the (C-NLS) equation, we utilize the improved modified, extended tanh-function method. By separating the Schrödinger equation into real and imaginary parts, we can obtain four coupled equations, which we then analyze using the generalized tanh method to extract exact solutions. This system of equations is essential for understanding the behavior of quantum systems and has various applications in quantum mechanics. We obtain an analytical solution and demonstrate numerical solutions using implicit finite difference. Studies have shown that this scheme is second-order in space and time, and the von Neumann stability analysis confirms its unconditional stability. We introduce the comparison between numerical and exact solutions.

1. Introduction

Nonlinear partial differential equations (PDEs) are crucial in describing various physical, biological, and engineering phenomena. The nonlinear Schrödinger equation (NLSE) is a pivotal partial differential equation in classical and quantum mechanics. The nonlinear Schrödinger equation (NLSE) is a partial differential equation that describes wave propagation dynamics in various physical systems, particularly in nonlinear optics and Bose–Einstein condensates. Therefore, the coupled nonlinear Schrödinger equations provide a robust framework for understanding and modeling the behavior of coherent optical pulses in nonlinear media. These equations are essential for investigating soliton dynamics, pulse compression, and mode interactions in various optical systems. The future behavior of these problems is well-known from the exact and numerical solutions of the corresponding PDEs. In other words, to better understand the long-term behavior of nonlinear phenomena, exploring their exact solutions is consequential. Therefore, developing fundamental and systematic methods for deriving analytical solutions to PDEs has become a popular and fascinating subject for most scholars. Among these techniques, we propose the Kudryashov approach [1,2], the improved $\mathcal{Q}$-expansion strategy [3,4], the $\left(\frac{G^{\prime }}{G}\right)$-expansion method [5], and the Jacobi elliptic expansion [6,7]. These methods are handy for transforming a given PDE into a more straightforward ordinary differential equation (ODE) that can be more easily solved. Furthermore, numerical methods are often required for solving nonlinear PDEs when analytical solutions are not feasible. Some commonly used numerical methods to solve nonlinear PDEs include the finite difference method (FDM) [8,9], the compact finite-difference method [10], the Galerkin finite element method (GFEM) [11], and the adaptive mesh refinement method, which is a powerful approach for efficiently solving PDEs while maintaining high accuracy, making it particularly valuable for simulations involving complex and dynamic phenomena [12,13,14,15,16], and more others [17,18,19]. The chosen method depends on the problem’s nature, the desired accuracy, computational resources, and available software libraries. Additionally, it is essential to consider stability, convergence, and efficiency when selecting a numerical method for solving nonlinear PDEs. Some other analytical and numerical solution methodologies can be found in the literature [20,21,22,23,24,25]. Several studies have examined the coupled nonlinear Schrödinger equation for its numerical and exact solutions. For example, references [26,27,28] introduced a finite difference method for a numerical simulation of the coupled nonlinear Schrodinger equation. Jianqiang Gu et al. [29] obtained analytical solutions via the collocation approach, the generalized Kudryashov method, the modified Kudryashov method, and the exponential rational function method. In [30], the extended, modified auxiliary equation mapping method establishes several solutions for the coupled nonlinear Schrödinger equations. Previous studies on the coupled nonlinear Schrödinger equations have mostly focused on finding the analytical or numerical solutions with the same type of equations in the system. Therefore, our study developed some traveling wave and numerical solutions of the system for two different types of equations. Consider the second-order coupled nonlinear Schrödinger (C-NLS) equation:

$$\begin{array}{cc}\hfill & \iota {\psi }_t+\iota \beta {\psi }_x-\alpha {\psi }_{xx}+\left({\left|\psi \right|}^2+\gamma {\left|\varphi \right|}^2\right)\psi =0,\hfill \\ & \iota {\varphi }_t+\alpha {\varphi }_{xx}+\left({\left|\varphi \right|}^2+\gamma {\left|\psi \right|}^2\right)\varphi =0.\hfill \end{array}$$

(1)

where $\psi (x,t)$ and $\varphi (x,t)$ are complex unknown functions, $\beta ,\alpha$ describes the dispersion in the optic fiber, $\gamma$ is the self-phase modulation parameter.

The initial conditions are as follows:

$$\psi (x,0)={\mathcal{F}}_1\left(x\right),\varphi (x,0)={\mathcal{F}}_2\left(x\right),$$

(2)

the following are the boundary conditions:

$${\psi }_x(x,t)=0,{\varphi }_x(x,t)=0,\mathrm{at}x=x_L,x_R.$$

(3)

We write the complex functions $\psi$ and $\varphi$ as a sum of their real and imaginary parts.

$$\psi (x,t)=u_1(x,t)+\iota u_2(x,t),\varphi (x,t)=v_1(x,t)+\iota v_2(x,t),$$

(4)

where $u_j,v_j$, $j=\left\{1,2\right\}$ are real-valued functions. By superseding Equation (4) into system (1), we obtain the following system

$$\begin{array}{cc}\hfill & \frac{\partial u_1}{\partial t}+\beta \frac{\partial u_1}{\partial x}-\alpha \frac{{\partial }^2{u}_2}{\partial x^2}+\left((u_1^2+u_2^2)+\gamma (v_1^2+v_2^2)\right)u_2=0,\hfill \\ & \frac{\partial u_2}{\partial t}+\beta \frac{\partial u_1}{\partial x}+\alpha \frac{{\partial }^2{u}_1}{\partial x^2}-\left((u_1^2+u_2^2)+\gamma (v_1^2+v_2^2)\right)u_1=0,\hfill \\ & \frac{\partial v_1}{\partial t}+\alpha \frac{{\partial }^2{v}_2}{\partial x^2}+\left((v_1^2+v_2^2)+\gamma (u_1^2+u_2^2)\right)v_2=0,\hfill \\ & \frac{\partial v_2}{\partial t}-\alpha \frac{{\partial }^2{v}_1}{\partial x^2}-\left((v_1^2+v_2^2)+\gamma (u_1^2+u_2^2)\right)v_1=0.\hfill \end{array}$$

(5)

Although numerous researchers work analytically to find the traveling wave solutions of the NLS equation Equation (1), only a few scientists investigate the analytical and numerical solutions of this problem after converting it into a system of real and imaginary parts. The generalized tanh method is a powerful technique for searching traveling waves from one-dimensional nonlinear waves and evolution equations to extract exact solutions of system (5). Hence, the numerical approach is employed to generate accurate and dependable results. One of the greatest ways to ensure accurate solutions is to check the conformity of exact and numerical solutions. Even though some experts only find exact solutions, in this study, we compare the exact and numerical solutions to ensure that the solutions are perfectly accurate and correct.

This paper is structured as follows: Section 2 describes the improved modified extended tanh-function method utilized to express the exact traveling wave solutions, and then apply it to extract solutions of system (1), while in Section 3, we obtain the analytical solution to (5) using the generalized tanh method. Section 4 is dedicated to solving system (5) by employing the implicit finite difference. Section 5 introduces the results and discussion. Section 6 highlights the most critical results discovered in this article.

2. Improved Modified Extended Tanh-Function Technique

This section established the improved modified extended tanh-function method and soliton solutions of nonlinear evolution equations (NLEEs). We assume that the following nonlinear evolution equations (NLEEs)

$${\Omega }_1\left({\psi }_t,{\psi }_x,{\psi }_{xx},{\varphi }_t,{\varphi }_x,{\varphi }_{xx},\dots \right)=0,$$

(6)

where $\psi =\psi (x,t),\varphi =\varphi (x,t)$.

The transformation of waves can be represented as

$$\psi (x,t)=u\left(\zeta \right)e^{\iota (x-\sigma t+ϵ)},\varphi (x,t)=v\left(\zeta \right)e^{\iota (x-\sigma t+ϵ)},\zeta =x-wt.$$

(7)

By substituting Equation (7) into Equation (6), we obtain the following ordinary differential equations (ODEs)

$${\Omega }_2\left(u^{\prime },u^{″},v^{\prime },v^{″},\dots \right)=0.$$

(8)

The traveling wave solution of Equation (8), according to the modified extended tanh-function technique [31], takes the following form

$$u\left(\zeta \right)=\sum _{j=0}^{N_1}{a}_j\Theta {\left(\zeta \right)}^j+\sum _{j=1}^{N_1}{\widehat{a}}_j\Theta {\left(\zeta \right)}^{-j},v\left(\zeta \right)=\sum _{j=0}^{N_2}{b}_j\Theta {\left(\zeta \right)}^j+\sum _{j=1}^{N_2}{\widehat{b}}_j\Theta {\left(\zeta \right)}^{-j}.$$

(9)

where $a_j,{\widehat{a}}_j,b_j$, and ${\widehat{b}}_j$ are evaluated later. The values of $N_1$ and $N_2$ can be obtained by finding the homogeneous balance of the nonlinear term and the highest derivative. The function $\Theta \left(\zeta \right)$ solves the following Riccati differential equation:

$${\Theta }^{\prime }=\lambda +{\Theta }^2.$$

(10)

Equation (10) has the general solutions [32], where $\lambda \ge 0$ or $\lambda <0$. The method is improved by generalizing Riccati Equation (10) to the following general ODE [33]:

$${\Theta }^{\prime }\left(\zeta \right)=\widehat{z}\sqrt{c_0+c_1\Theta +c_2{\Theta }^2+c_3{\Theta }^3+c_4{\Theta }^4}.$$

(11)

We have a set of constants, denoted by $c_l$, where l belongs to set $\{0,1,2,3,4\}$ and subject to certain restrictions. Additionally, $\widehat{z}$ is equal to either 1 or $-1$. We need to determine the values of $N_1$ and $N_2$ for use in Equation (9). After that, we can combine Equations (9) and (11) with Equation (8) to obtain algebraic equations. We obtain a system of equations by collecting the coefficients for the same power of $\Theta \left(\zeta \right)$. Solving this system gives us the values of $a_j$, ${\widehat{a}}_j$, $b_j$, and ${\widehat{b}}_j$.

Application of the Method

New traveling wave solutions of system (1) are obtained by substituting Equation (7) into (1). The coefficients are then collected and equated to zero for both imaginary and real parts, resulting in some relations:

$$\begin{array}{cc}\hfill & \left(-2\alpha +\beta -w\right)u^{\prime }=0,\hfill \\ & -\alpha u^{″}+(\sigma -\beta +\alpha )u++\gamma u{v}^2+u^3=0,\hfill \\ & \left(2\alpha -w\right)v^{\prime }=0,\hfill \\ & \alpha v^{″}+(\sigma -\alpha )v++\gamma v{u}^2+v^3=0.\hfill \end{array}$$

(12)

We obtain w and $\beta$ as follows:

$$w=2\alpha ,\beta =4\alpha .$$

(13)

When the conditions Equation (13) are applied, then Equation (12) becomes

$$\begin{array}{cc}\hfill & -\alpha u^{″}+(\sigma -3\alpha )u+\gamma u{v}^2+u^3=0,\hfill \\ & \alpha v^{″}+(\sigma -\alpha )v+\gamma v{u}^2+v^3=0.\hfill \end{array}$$

(14)

By using the balancing procedure between $u^{″}$ and $u^3$ in the first equation, and between $v^{″}$ and $v^3$ in the second equation, we obtain $N_1=N_2=1$. Equation (9) takes the form

$$u=a_0+a_1\Theta +{\widehat{a}}_1{\Theta }^{-1},v=b_0+b_1\Theta +{\widehat{b}}_1{\Theta }^{-1}.$$

(15)

To solve $a_j$, $b_j$, ${\widehat{a}}_1$, ${\widehat{b}}_1$, $\gamma$, $\sigma$, we substitute (11) and (15) into (14) and collect terms with the same order of $\Theta$. This converts the left-hand side of (14) into a polynomial in $\Theta$. Equating each coefficient of the resulting polynomial to zero gives us a set of algebraic equations. We can solve these equations using Mathematica and obtain the following values, where $j=\left\{0,1\right\}$:

$\mathrm{If}{c}_1=c_3=0,c_0=\frac{c_2^2}{4c_4}$, then

$$\begin{array}{cc}& a_1=\frac{1}{2}\sqrt{\pm \frac{s}{c_0}+4\alpha c_4{\widehat{z}}^2},b_1=\pm \frac{1}{2}\sqrt{\pm \frac{s}{c_0}-4\alpha c_4{\widehat{z}}^2},{\widehat{a}}_1=\pm \frac{s\sqrt{\pm \frac{s}{c_0}+4\alpha c_4{\widehat{z}}^2}}{2\alpha c_4\left(c_2{\widehat{z}}^2+1\right)},\hfill \\ & {\widehat{b}}_1=-\frac{s\sqrt{\frac{s-4\alpha c_0{c}_4{\widehat{z}}^2}{c_0}}}{2\alpha c_4\left(c_2{\widehat{z}}^2+1\right)},\sigma =\frac{2\left(\alpha +\alpha c_2{\widehat{z}}^2\mp s{\widehat{z}}^2\right)}{c_2{\widehat{z}}^2+1},s=\sqrt{c_0{c}_4\left(\alpha +\alpha c_2{\widehat{z}}^2\right){}^2},\hfill \\ & \kappa =\sqrt{{\alpha }^2{c}_2^2\left(c_2+1\right){}^2},\gamma =-1,a_0=b_0=0.\hfill \end{array}$$

Therefore, the hyperbolic function solutions for system (1) are as follows:

$$\begin{array}{cc}\hfill {\psi }_1(x,t)=& \frac{\left(\sqrt{-\frac{c_2}{c_4}}\sqrt{\frac{c_4\left(2\alpha c_2^2+\kappa \right)}{c_2^2}}\tanh \left(\frac{\sqrt{-c_2}(x-2\alpha t)}{\sqrt{2}}\right)\left(\alpha c_2^2+\alpha c_2-\kappa {\coth }^2\left(\frac{\sqrt{-c_2}(x-2\alpha t)}{\sqrt{2}}\right)\right)\right)}{2\alpha c_2\left(c_2+1\right)}\hfill \\ & \times \exp \left(\frac{\iota \left(c_2(-2\alpha t+x+ϵ)-2\alpha t+\kappa t+x+ϵ\right)}{c_2+1}\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\varphi }_1(x,t)=& \frac{\sqrt{-\frac{c_2}{c_4}}\sqrt{\frac{c_4\left(\kappa -2\alpha c_2^2\right)}{c_2^2}}\tanh \left(\frac{\sqrt{-c_2}(x-2\alpha t)}{\sqrt{2}}\right)\left(\alpha c_2^2+\alpha c_2+\kappa {\coth }^2\left(\frac{\sqrt{-c_2}(x-2\alpha t)}{\sqrt{2}}\right)\right)}{2\alpha c_2\left(c_2+1\right)}\hfill \\ & \times \exp \left(\frac{\iota \left(c_2(-2\alpha t+x+ϵ)-2\alpha t+\kappa t+x+ϵ\right)}{c_2+1}\right).\hfill \end{array}$$

(17)

Then, the trigonometric function solutions of system (1) can be written as

$$\begin{array}{cc}\hfill {\psi }_2(x,t)=& \frac{\sqrt{\frac{c_2}{c_4}}\sqrt{\frac{c_4\left(2\alpha c_2^2+\kappa \right)}{c_2^2}}\tan \left(\frac{\sqrt{c_2}(x-2\alpha t)}{\sqrt{2}}\right)\left(\alpha c_2^2+\alpha c_2\pm \kappa {\cot }^2\left(\frac{\sqrt{c_2}(x-2\alpha t)}{\sqrt{2}}\right)\right)}{2\alpha c_2\left(c_2+1\right)}\hfill \\ & \times \exp \left(\frac{\iota \left(c_2(-2\alpha t+x+ϵ)-2\alpha t+\kappa t+x+ϵ\right)}{c_2+1}\right),\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\varphi }_2(x,t)=& \pm \frac{\sqrt{\frac{c_2}{c_4}}\sqrt{\frac{c_4\left(\kappa -2\alpha c_2^2\right)}{c_2^2}}\tan \left(\frac{\sqrt{c_2}(x-2\alpha t)}{\sqrt{2}}\right)\left(\alpha c_2^2+\alpha c_2\mp \kappa {\cot }^2\left(\frac{\sqrt{c_2}(x-2\alpha t)}{\sqrt{2}}\right)\right)}{2\alpha c_2\left(c_2+1\right)}\hfill \\ & \times \exp \left(\frac{\iota \left(c_2(-2\alpha t+x+ϵ)-2\alpha t+\kappa t+x+ϵ\right)}{c_2+1}\right).\hfill \end{array}$$

(19)

Figure 1 and Figure 2 show the real and imaginary parts of ${\psi }_1(x,t)$ and ${\varphi }_1(x,t)$, respectively. These surfaces clarify the physical meaning of the discussed equation. We have $w,\beta$, and we take it, depending on the condition of Equation (13).

3. Generalized Tanh Method

We simply describe the fundamental steps of the generalized tanh method [32]. The generalized tanh method replaces the tanh function in the tanh method with a Riccati equation solution. Consider the general form of the nonlinear partial differential equations (NPDEs) as follows:

$${\Omega }_1\left(u_{1,t},u_{1,x},u_{1,xx},u_{2,t},u_{2,x},u_{2,xx},v_{1,t},v_{1,x},v_{1,xx},v_{2,t},v_{2,x},v_{2,xx},\dots \right)=0,$$

(20)

where $u_j=u_j(x,t),v_j=v_j(x,t),j=\left\{1,2\right\}$. Then we employ the following wave transformation

$$\begin{array}{cc}\hfill & u_1=U_1\left(\zeta \right),u_2=U_2\left(\zeta \right),v_1=V_1\left(\zeta \right),v_2=V_2\left(\zeta \right),\zeta =x-wt.\hfill \end{array}$$

(21)

By substituting Equation (21) into Equation (20), we obtain the following ordinary differential equations (ODEs)

$${\Omega }_2\left(U_1^{\prime },U_1^{″},U_2^{\prime },U_2^{″},V_1^{\prime },V_1^{″},V_2^{\prime },V_2^{″},\dots \right)=0.$$

(22)

Then we propose the following series expansion as a solution to Equation (22):

$$\begin{array}{cc}\hfill & U_1=a_0+\sum _{j=1}^{N_1}{a}_j{\Phi }^j,U_2=c_0+\sum _{j=1}^{N_2}{c}_j{\Phi }^j,\hfill \\ & V_1=b_0+\sum _{j=1}^{N_3}{b}_j{\Phi }^j,V_2=d_0+\sum _{j=1}^{N_4}{d}_j{\Phi }^j.\hfill \end{array}$$

(23)

where $\Phi =\Phi \left(\zeta \right)$ satisfies the following ordinary differential equation

$${\Phi }^{\prime }=\varpi +{\Phi }^2.$$

(24)

The positive integers $N_l,l=\left\{1,2,3,4\right\}$ can be determined by balancing the highest derivative term with nonlinear terms in (22). Substituting (23) and (24) into (22) and then setting zero all coefficients for the same order of $\Phi$, we can obtain a system of algebraic equations, from which the constants $a_0,a_j,b_0$, and $b_j$ are obtained explicitly. The Riccati equation admits various kinds of solutions, presented in [32]. By using the solution to Equation (24), we obtain:

$$\Phi =\left\{\begin{array}{c}-\sqrt{-\varpi }\tanh \left(\sqrt{-\varpi }\zeta \right)\hfill \\ -\sqrt{-\varpi }\coth \left(\sqrt{-\varpi }\zeta \right)\hfill \end{array}\right.\mathrm{for}\varpi <0.$$

(25)

Application of the Method

To find the solitary wave solutions of system (5), we use the wave transformation (21), the system (5) can be converted to the following ODEs:

$$\begin{array}{cc}\hfill & -w{U}_1^{\prime }+\beta U_1^{\prime }-\alpha U_2^{″}+\left(\left(U_1^2+U_2^2\right)+\gamma \left(V_1^2+V_2^2\right)\right)U_2=0,\hfill \\ & -w{U}_2^{\prime }+\beta U_2^{\prime }+\alpha U_1^{″}+\left(\left(U_1^2+U_2^2\right)+\gamma \left(V_1^2+V_2^2\right)\right)U_1=0,\hfill \\ & -w{V}_1^{\prime }+\alpha V_2^{″}+\left(\left(V_1^2+V_2^2\right)+\gamma \left(U_1^2+U_2^2\right)\right)V_2=0,\hfill \\ & -w{V}_2^{\prime }-\alpha V_1^{″}+\left(\left(V_1^2+V_2^2\right)+\gamma \left(U_1^2+U_2^2\right)\right)V_1=0.\hfill \end{array}$$

(26)

By balancing the highest derivative and nonlinear terms for each equation in system (26), we obtain $N_l=1$, where $l=\left\{1,2,3,4\right\}$, so according to Equation (23), we assume that

$$\begin{array}{cc}\hfill & U_1=a_0+a_1\Phi ,U_2=c_0+c_1\Phi ,\hfill \\ & V_1=b_0+b_1\Phi ,V_2=d_0+d_1\Phi .\hfill \end{array}$$

(27)

Substituting Equation (27) and the necessary derivatives into system (26) using Equation (25), applying trigonometric identities and collecting the coefficients of the same power of $\Phi$ that contain independent combinations to zero, we obtain the following cases and solutions

$$\begin{array}{cc}\hfill \mathbf{Family}\mathbf{1}:& a_0=\sqrt{d_0^2-2\alpha },b_1=\mp \sqrt{d_0^2-2\alpha },c_0=\mp \sqrt{d_0^2-2\alpha },\hfill \\ & d_1=\mp \sqrt{d_0^2-2\alpha },w=\pm \frac{2\alpha d_0}{\sqrt{d_0^2-2\alpha }},\beta =\pm \frac{4\alpha \left(d_0^2-\alpha \right)}{d_0\sqrt{d_0^2-2\alpha }},\hfill \\ & \gamma =-1,\varpi =-1,a_1=\pm d_0,c_1=d_0,b_0=-d_0.\hfill \end{array}$$

Thus, the new exact solutions of the (C-NLS) can be identified as:

$$\begin{array}{cc}\hfill & u_1(x,t)=\sqrt{d_0^2-2\alpha }\mp d_0\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right),\hfill \\ & u_2(x,t)=\mp \sqrt{d_0^2-2\alpha }-d_0\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right),\hfill \\ & v_1(x,t)=\pm \sqrt{d_0^2-2\alpha }\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right)-d_0,\hfill \\ & v_2(x,t)=\pm \sqrt{d_0^2-2\alpha }\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right)+d_0.\hfill \end{array}$$

(28)

Thus, by substituting Equation (28) into Equation (4), we obtain the exact solution to system (1):

$$\begin{array}{cc}\hfill & {\psi }_3(x,t)=\sqrt{d_0^2-2\alpha }\mp d_0\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right)+\iota \left(\mp \sqrt{d_0^2-2\alpha }-d_0\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right)\right),\hfill \\ \hfill & {\varphi }_3(x,t)=-d_0\pm \sqrt{d_0^2-2\alpha }\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right)+\iota \left(d_0\pm \sqrt{d_0^2-2\alpha }\tanh \left(x-\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}\right)\right).\hfill \end{array}$$

(29)

The traveling wave solutions of (C-NLS):

$$\begin{array}{cc}\hfill & u_1(x,t)=\sqrt{d_0^2-2\alpha }\pm d_0\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right),\hfill \\ & u_2(x,t)=\pm \sqrt{d_0^2-2\alpha }-d_0\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right),\hfill \\ & v_1(x,t)=\mp d_0\mp \sqrt{d_0^2-2\alpha }\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right),\hfill \\ & v_2(x,t)=d_0\mp \sqrt{d_0^2-2\alpha }\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right).\hfill \end{array}$$

(30)

Thus, by substituting Equation (30) into Equation (4), we obtain the exact solution to system (1):

$$\begin{array}{cc}\hfill & {\psi }_4(x,t)=\sqrt{d_0^2-2\alpha }\pm d_0\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right)+\iota \left(\pm \sqrt{d_0^2-2\alpha }-d_0\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right)\right),\hfill \\ & {\varphi }_4(x,t)=\mp d_0\mp \sqrt{d_0^2-2\alpha }\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right)+\iota \left(d_0\mp \sqrt{d_0^2-2\alpha }\coth \left(\frac{2\alpha d_{0}t}{\sqrt{d_0^2-2\alpha }}+x\right)\right).\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \mathbf{Family}\mathbf{2}:& a_0=2\sqrt{\alpha },a_1=\pm 2\sqrt{\alpha },b_0=\sqrt{2}\sqrt{\alpha },b_1=\pm \sqrt{2}\sqrt{\alpha },\hfill \\ & c_0=2\sqrt{\alpha },c_1=\mp 2\sqrt{\alpha },d_0=\sqrt{2}\sqrt{\alpha },d_1=\mp \sqrt{2}\sqrt{\alpha },\hfill \\ & w=\pm 2\alpha ,\beta =\pm 4\alpha ,\varpi =-1,\gamma =-1.\hfill \end{array}$$

Based on Family 2, we can identify the new exact solutions of (C-NLS):

$$\begin{array}{cc}\hfill & u_1(x,t)=\mp 2\sqrt{\alpha }\left(\tanh \left(x-2\alpha t\right)\mp 1\right),u_2(x,t)=\pm 2\sqrt{\alpha }\left(\tanh \left(x-2\alpha t\right)\pm 1\right),\hfill \\ & v_1(x,t)=\mp \sqrt{2}\sqrt{\alpha }\left(\tanh \left(x-2\alpha t\right)\mp 1\right),v_2(x,t)=\pm \sqrt{2}\sqrt{\alpha }\left(\tanh \left(x-2\alpha t\right)\pm 1\right).\hfill \end{array}$$

(32)

Thus, by substituting Equation (32) into Equation (4), we obtain the exact solution to system (1):

$$\begin{array}{cc}\hfill & {\psi }_5(x,t)=\mp 2\sqrt{\alpha }\left(\tanh \left(x-2\alpha t\right)\mp 1\right))\pm 2\sqrt{\alpha }\iota \left(\tanh \left(x-2\alpha t\right)\pm 1\right),\hfill \\ \hfill & {\varphi }_5(x,t)=\mp \sqrt{2\alpha }\left(\tanh \left(x-2\alpha t\right)\mp 1\right)\pm \sqrt{2\alpha }\iota \left(\tanh \left(x-2\alpha t\right)\pm 1\right).\hfill \end{array}$$

(33)

The traveling wave solutions of (C-NLS):

$$\begin{array}{cc}\hfill & u_1(x,t)=\mp 2\sqrt{\alpha }(\coth (x-2\alpha t)\mp 1),u_2(x,t)=\pm 2\sqrt{\alpha }(\coth (x-2\alpha t)\pm 1),\hfill \\ & v_1(x,t)=\mp \sqrt{2\alpha }(\coth (x-2\alpha t)\mp 1),v_2(x,t)=\pm \sqrt{2\alpha }(\coth (x-2\alpha t)\pm 1).\hfill \end{array}$$

(34)

Thus, by substituting Equation (34) into Equation (4), we obtain the exact solution to system (1):

$$\begin{array}{cc}\hfill & {\psi }_6(x,t)=\mp 2\sqrt{\alpha }(\coth (x-2\alpha t)\mp 1)\pm 2\sqrt{\alpha }\iota (\coth (x-2\alpha t)\pm 1),\hfill \\ \hfill & {\varphi }_6(x,t)=\mp \sqrt{2\alpha }(\coth (x-2\alpha t)\mp 1)\pm \sqrt{2\alpha }\iota (\coth (x-2\alpha t)\pm 1).\hfill \end{array}$$

(35)

Figure 3 and Figure 4 present the 3D surfaces of Equation (28) with specific parameter values. Figure 5 presents the behavior for the analytical solutions in Equation (28) when we change the values of $\alpha$, while the other parameters take fixed values $d_0=-2,N_x=1000$ with $t=1,x=-10\to 10$. As a result, we can see that changing the parameter $\alpha$ has a different effect between $u_1(x,t)$ and $u_2(x,t)$, while $v_1(x,t),v_2(x,t)$ has the same effect.

4. Numerical Solution

In this section, we employ the implicit finite-difference method to extract the numerical results of system (5) on a physical domain $[x_L,x_R]$. The domain is divided into $N_x$ subintervals $[x_m,x_{m+1}]$, such that

$$x_m=(m-1){\Delta }_x,\forall x_m\in [x_L,x_R],m=1,2,\dots ,N_x+1,$$

where ${\Delta }_x=(x_R-x_L)/N_x$, ${\Delta }_x$ denotes a uniform width of each subinterval. The implicit finite difference method is a numerical technique used to solve partial differential equations (PDEs) by discretizing the PDE in time and space. Unlike explicit methods, where the future solution is computed using information only from the current time step, implicit methods consider information from both the current and future time steps. This makes them unconditionally stable and suitable for solving stiff PDEs, where explicit methods might be impractical due to restrictive stability conditions.

The system presented in Equation (5) can be expressed using matrices and vectors, if we assume, ${\mathbf{u}}^{\mathrm{T}}=\left[\begin{array}{cc}{u}_1& u_2\end{array}\right],{\mathbf{v}}^{\mathrm{T}}=\left[\begin{array}{cc}{v}_1& v_2\end{array}\right]$, then

$$\frac{\partial \Gamma }{\partial t}+\beta A\frac{\partial \Gamma }{\partial x}-\alpha B\frac{{\partial }^2\Gamma }{\partial x^2}+\mathrm{G}\left(\Gamma \right)\Gamma =0,$$

(36)

where, ${\Gamma }^{\mathrm{T}}=\left[\begin{array}{cc}\mathbf{u}& \mathbf{v}\end{array}\right],A=\left[\begin{array}{cc}\mathrm{I}& z\\ z& z\end{array}\right],k_1=\left({\mathbf{u}}^{\mathrm{T}}\mathbf{u}+\gamma {\mathbf{v}}^{\mathrm{T}}\mathbf{v}\right),k_2=\left({\mathbf{v}}^{\mathrm{T}}\mathbf{v}+\gamma {\mathbf{u}}^{\mathrm{T}}\mathbf{u}\right)$, while $\mathrm{I},z\in {\mathbb{R}}^{2\times 2}$ identity matrix and zero matrix, respectively.

$$B=\left[\begin{array}{cc}\begin{array}{cc}0& 1\\ -1& 0\end{array}& z\\ z& \begin{array}{cc}0& -1\\ 1& 0\end{array}\end{array}\right],\mathrm{G}\left(\Gamma \right)=\left[\begin{array}{cc}\begin{array}{cc}0& k_1\\ -k_1& 0\end{array}& z\\ z& \begin{array}{cc}0& k_2\\ -k_2& 0\end{array}\end{array}\right].$$

Then, we use the approximate solutions ${\Gamma }_m$ to the analytical solution $\Gamma (x_m,t)=\Gamma$. To approximate the first and second derivatives in Equation (36), we use the following central difference formulas:

$$\begin{array}{cc}\hfill & \frac{\partial \Gamma }{\partial x}=\frac{1}{2{\Delta }_x}\left({\Gamma }_{m+1}-{\Gamma }_{m-1}\right)=\frac{1}{2{\Delta }_x}{\delta }_x{\Gamma }_m,\hfill \\ & \frac{{\partial }^2\Gamma }{\partial x^2}=\frac{1}{{\Delta }_x^2}\left({\Gamma }_{m+1}-2{\Gamma }_m+{\Gamma }_{m-1}\right)=\frac{1}{{\Delta }_x^2}{\delta }_x^2{\Gamma }_m.\hfill \end{array}$$

(37)

When we insert Equation (37) into Equation (36), we obtain the semi-discrete system

$${\Gamma }_t{|}_m+\frac{\beta }{2{\Delta }_x}A{\delta }_x{\Gamma }_m-\frac{\alpha }{{\Delta }_x^2}B{\delta }_x^2{\Gamma }_m+\mathrm{G}\left({\Gamma }_m\right){\Gamma }_m=0.$$

(38)

Equation (38) can be expressed in a vector form, as follows:

$${\Gamma }_t{|}_m+\mathrm{F}\left({\Gamma }_m\right)=0,$$

(39)

where $\mathrm{F}\left({\Gamma }_m\right)=\frac{\beta }{2{\Delta }_x}A{\delta }_x{\Gamma }_m-\frac{\alpha }{{\Delta }_x^2}B{\delta }_x^2{\Gamma }_m+\mathrm{G}\left({\Gamma }_m\right){\Gamma }_m$.

We solve system (39) by applying the implicit midpoint rule. Here, ${\Gamma }_m^n$ represents the discrete approximation of the exact solution $\Gamma (x_m,t_n)$. The semi-discretization of (39) is obtained as follows:

$${\Gamma }_t{|}_m^n+\mathrm{F}\left(\frac{{\Gamma }_m^{n+1}+{\Gamma }_m^n}{2}\right)=0.$$

(40)

The above system was solved using an ODE solver in FORTRAN, denoted as the DDASPK solver [34]. This solver uses implicit differentiation operators to approximate time derivatives. The obtained numerical results are acceptable. The Taylor expansion is utilized to determine the accuracy of the numerical scheme. When plugging the expansion and simplifying the equation, the accuracy is $\mathcal{O}({\Delta }_t^2,{\Delta }_x^2)$. This means the scheme has a second-order accuracy in both space and time. Then, we analyze the stability of the numerical method using the von Neumann stability analysis. This analysis only applies to linear difference schemes, so we examine the linearized form of Equation (40).

$${\Gamma }_m^{n+1}-{\Gamma }_m^n+{\Delta }_t\mathcal{Q}\left(\frac{{\Gamma }_m^{n+1}+{\Gamma }_m^n}{2}\right),$$

(41)

where

$$\mathcal{Q}\left({\Gamma }_m^n\right)=\frac{\beta }{2{\Delta }_x}A{\delta }_x{\Gamma }_m^n-\frac{\alpha }{{\Delta }_x^2}B{\delta }_x^2{\Gamma }_m^n+\tilde{A}{\Gamma }_m^n,$$

$\tilde{A}=\left[\begin{array}{cc}\begin{array}{cc}0& {\tilde{k}}_1\\ -{\tilde{k}}_1& 0\end{array}& z\\ z& \begin{array}{cc}0& {\tilde{k}}_2\\ -{\tilde{k}}_2& 0\end{array}\end{array}\right],\mathrm{and}{\tilde{k}}_1=max\left\{k_1\right\},{\tilde{k}}_2=max\left\{k_2\right\}.$ To analyze the stability of the difference scheme, we assume

$${\Gamma }_m^n={\mathbf{M}}^n{e}^{\iota \omega m{\Delta }_x}.$$

(42)

When using von Neumann, the following criterion is required for stability.

$$\underset{j}{max}\left|{\mu }_j\right|\le 1,\mathrm{where}j=1,2,3,4.$$

After substituting Equation (42) into Equation (41), we obtain the amplification matrix as a result of a calculation, which can be expressed as a matrix equation:

$$\mathbf{M}\left(\mathbf{I}+\iota {\vartheta }_{1}A+{\vartheta }_{2}B+\frac{{\Delta }_t}{2}\tilde{A}\right)=\left(\mathbf{I}-\iota {\vartheta }_{1}A-{\vartheta }_{2}B-\frac{{\Delta }_t}{2}\tilde{A}\right),$$

(43)

where ${\vartheta }_1=\frac{\beta {\Delta }_t}{2{\Delta }_x},{\vartheta }_2=\frac{2\alpha {\Delta }_t}{{\Delta }_x^2}{\sin }^2\left(\frac{\omega {\Delta }_x}{2}\right).$

The matrix M can be given explicitly as

$$\mathbf{M}={\left(\mathbf{I}+\iota {\vartheta }_{1}A+{\vartheta }_{2}B+\frac{{\Delta }_t}{2}\tilde{A}\right)}^{-1}\left(\mathbf{I}-\iota {\vartheta }_{1}A-{\vartheta }_{2}B-\frac{{\Delta }_t}{2}\tilde{A}\right).$$

(44)

We calculated the eigenvalues of the matrix M, and they take the following forms:

$$\begin{array}{cc}\hfill & {\mu }_1=\frac{-\left({\vartheta }_4-{\vartheta }_2\right)+\iota }{\left({\vartheta }_4-{\vartheta }_2\right)+\iota },{\mu }_2=\frac{\left({\vartheta }_4-{\vartheta }_2\right)+\iota }{-\left({\vartheta }_4-{\vartheta }_2\right)+\iota },\hfill \\ & {\mu }_3=-\frac{\left({\vartheta }_2+{\vartheta }_3\right)+{\vartheta }_1+\iota }{\left({\vartheta }_2+{\vartheta }_3\right)+{\vartheta }_1-\iota },{\mu }_4=\frac{-\left({\vartheta }_2+{\vartheta }_3\right)+{\vartheta }_1+\iota }{\left({\vartheta }_2+{\vartheta }_3\right)-{\vartheta }_1+\iota },\hfill \end{array}$$

(45)

where

$${\vartheta }_3=\frac{{\Delta }_t}{2}{\tilde{k}}_1,{\vartheta }_4=\frac{{\Delta }_t}{2}{\tilde{k}}_2.$$

The modulus of these eigenvalues is equal to 1

$$|{\mu }_j|=1,j=1,2,3,4.$$

(46)

The stability condition of the von Neumann analysis is satisfied. The scheme is unconditionally stable.

5. Results and Discussion

Using the improved modified extended tanh-function method, we obtained some exact closed-form solutions for the second-order coupled nonlinear Schrödinger (C-NLS) Equation (1). Figure 1 and Figure 2 show the exact solution ${\psi }_1(x,t),{\varphi }_1(x,t)$ of system (1), where Figure 1a,b present real and imaginary parts, respectively. While Figure 2a presents real parts and Figure 2b presents imaginary parts of ${\varphi }_1(x,t)$. The parameter values are $\alpha =-4,ϵ=1,$ $c_2=-2,c_4=1, z=1,N_x=1000$ with $t=0\to 10$ and $x=(-5,3)$. After converting the coupled (NLS) equation to the system of real and imaginary Equations (5), we investigate the analytical solutions of system (5) using the generalized tanh method. Figure 3 and Figure 4 present the 3D surfaces of Equation (28) using the parameter values $\alpha =-1/8,$ $d_0=-1/2,N_x=1000$ with $t=0\to 10,x=-10\to 10$. Figure 5 presents the behavior for the analytical solutions in Equation (28) when we change the values of $\alpha$, while the other parameters take fixed values $d_0=-2,N_x=1000$ with $t=1,x=-10\to 10$. As a result, we can see that changing the parameter $\alpha$ has a different effect between $u_1(x,t)$ and $u_2(x,t)$, while $v_1(x,t),v_2(x,t)$ has the same effect. We examine numerical solutions using the implicit finite difference method to convert the underlying problems into a system of ODEs while maintaining continuous time derivatives. Figure 6 and Figure 7 and

Table 1. At time $t=10$, this is the relative error with the $L_2$ norm and CPU.

${\mathit{N}}_{\mathit{x}}$	The Relative Error	CPU
$100$	$6.50\times {10}^{-1}$	$0.043\times {10}^3\mathrm{s}$
200	$1.30\times {10}^{-2}$	$0.16\times {10}^3$ s
400	$5.50\times {10}^{-3}$	$0.41\times {10}^3$ s
800	$2.20\times {10}^{-4}$	$0.91\times {10}^3$ s
1000	$4.20\times {10}^{-5}$	$2.81\times {10}^3$ s
2000	$1.32\times {10}^{-5}$	$5.50\times {10}^3$ s

present the most important aspects of the solutions, enabling a direct comparison between the traveling wave solutions obtained by using the proposed exact methods with $N_x=2000$ and the numerical solutions obtained by using the finite difference method with various mesh numbers in the x direction. The ability to make this comparison is facilitated by the data presented in Figure 6 and Figure 7, as well as in Table 1. The numerical outcomes are comparable to a significant degree as a consequence of this. When ${\Delta }_x$ scores reach zero, the mean error also reaches zero. The numerical methods are invariably stable when the parameter values are set to $\alpha =-0.125,d_0=-0.5.$ This method produces reliable and powerful results (Figure 8).

6. Conclusions

In this study, we achieve exact traveling wave solutions of the coupled nonlinear Schrödinger (C-NLS) Equation (1) using an improved modified extended tanh-function method. Furthermore, we convert the complex-valued function into its real and imaginary parts into a system of differential equations over a real field to solve the nonlinear Schrödinger equation in one dimension. Also, we investigate the exact and numerical solutions for system (5) using the generalized tanh method and implicit finite difference scheme, respectively. The numerical solutions approximately approach the exact solutions for small ${\Delta }_x$. More precisely, the $L_2$ error rapidly declines for smaller ${\Delta }_x$, as presented in Table 1 and Figure 8. The numerical scheme is found unconditionally stable via von Neumann stability analysis, exhibiting second-order accuracy in both time and space. We use FORTRAN 95 software by applying ODE solvers to extract approximate solutions to the problem and use Mathematica 13.2 software to investigate the exact solutions. Therefore, we use MATLAB R2023b software to create 3D graphs that accurately depict the solution based on specific parameter values. Also, in Figure 5, we plot 2D figures to illustrate the solution behavior when changing $\alpha$ and fixing others. Furthermore, the 3D figures illustrate that the exact and numerical solutions coincide and agree. The techniques utilized in this study can be applied to other NLPDEs in natural sciences.