In recent years the concept of soliton has been receiving considerable attention in optical communications. Since soliton is capable of propagating over long distances without change of shape and velocity, it has been found that the soliton propagating through optical fiber arrays is governed by a set of equations related to the coupled nonlinear Schrödinger equation   .
where , is the envelope or the amplitude of the jth wave packets. Equation (1) reduces to the standard nonlinear Schrödinger equation for , to Manakov integrable system for , and recently for this case the exact two soliton solution obtained and novel shape changing in elastic collision property has been brought out. The system for is of physical interest, in optical communication, and in biophysics, this system can be used to study the lunching and propagation of solitons along the three spines of an alpha-helix shape changing in protein    . In this work, we are going to derive a numerical solution for the three coupled nonlinear Schrödinger equations
with initial conditions
and the homogenous boundary conditions
The exact soliton solution of the 3-coupled nonlinear Schrödinger equation  , is given by
where are four arbitrary complex parameters. Further gives the amplitude of the jth mode and the soliton velocity.
The proposed system is of physical interest, in optical communication, and in biophysics. This system can be used to study the lunching and propagation of solitons along the three spines of an alpha-helix shape changing in protein    . In this work we are going to derive a numerical method of sixth order in space and second order in time for the three coupled nonlinear Schrödinger Equations (2)-(4).
Many numerical methods for solving the coupled nonlinear Schrödinger equation are derived in the last two decades. Finite difference and finite element methods are used to solve this system by Ismail       . A conservative compact finite difference schemes are given in  . Xing Lü studied the bright soliton collisions with shape change by intensity for the coupled Sasa-Satsuma system in the optical fiber communications in  and . To avoid complex computations, we assume
where are real functions, by separating the real and imaginary parts, and we write,
and we have assumed
where are real functions.
by substituting (8) into(2)-(4), the following system is obtained:
the system (2)-(4) can be written in a matrix-vector form as
Proposition 1: The three coupled nonlinear Schrödinger equations have the conserved quantities
To prove the first conserved quantity (17), we have
by multiplying (21) by and (22) by , and by adding the resulting equations to obtain
Integrating Equation (23) with respect to x from xL to xR and using the vanishing boundary conditions to obtain
and this is the proof of the first conserved quantity (17) The other two conserved quantities (18) and (19) can be proved in the same way.
The exact values of the conserved quantities using the exact soliton solution (7) are given by the following formula
The paper is organized as follows. In Section 2, we derived the high order compact finite difference scheme. The Fixed-Point scheme is derived in Section 3 to solve the block nonlinear penta-diagonal systems obtained in Section 2. In Section 4, we study the stability of our scheme. The numerical results of the derived method are reported in Section 6. Finally, we draw some conclusions in Section 7.
The scheme in (33)-(38) is of sixth order accuracy in space and second order in time, and it is unconditionally stable using von-Neumann stability analysis. A nonlinear block tridiagonal system must be solved at each time step. Fixed point method is used to do this job, and this will be discussed later.
2. High Order Compact Finite Difference Scheme
The compact finite difference is a numerical method to compute finite difference approximations. Such approximations tend to be more accurate for their stencil size (i.e. their compactness) and, for hyperbolic problems, have favorable dispersive error and dissipative error properties when compared to explicit schemes . In order to develop a numerical method for solving the system given in (2)-(4), the region will be covered with a rectangular mesh of points with coordinates,
where h and k are the space and time increments respectively. We denote the exact and numerical solution at the grid point by and , respectively. To evaluate the second derivatives at interior nodes, we assume that they can be obtained by solving the following penta-diagonal system 
Now, by Taylor Expansion, we can have the truncation error as the following
if we solve and , we get
so, the truncation error becomes
if then and which gives the explicit fourth-order scheme for the second derivative. Furthermore, when , the scheme becomes sixth order accurate, in this case and . By substituting these on formula (25) and after simplification, the space derivative of sixth order can be given implicitly as
Imposing the approximation given on the spatial direction, by using (9)-(14) into Equation (26), we get
The Crank-Nicolson discretization on the temporal direction of the 3-CNLS equation to obtain the numerical scheme
Equations (33)-(38) form a block pentadiagonal system as the following
The present method is of second order accuracy in time and sixth order in space, it is unconditionally stable, see Ismail . The resulting system is a block nonlinear penta-diagonal system which can be solved by fixed point method and this will be discussed next.
3. Fixed Point Method
Since the compact finite difference scheme (40) is nonlinear and implicit, an iterative method is needed to solve it. The fixed point for solving the resulting system can be given in a matrix vector form as follows  .
where the superscript s denotes the sth iterate for solving the nonlinear system of equations for each iteration. The system in (41) can be solved by Crout’s method, where we need only one LU factorization for the block-pentadiagonal matrix at the beginning of the calculation, and the solutions of lower and upper pentadiagonal block systems at each iteration are required only. The initial iterate can be chosen as . We apply the iterative schemes till the following condition
To study the stability of the scheme (2)-(4) we use the Von Neumann method, we do the following.
As we know Von Neumann method can be applied only for linear schemes, so we must study the linear version of (2)-(4) by freezing the nonlinear terms by assuming
The linear version of imposing the approximation given on the spatial direction (26) and the Crank-Nicolson discretization on the temporal direction of the Equations (2)-(4) can be displayed as follows:
where and .
By using (43) we can deduce the following relations
By substituting (43) and (44) into (42), we can get
We can write equations (45) as
So, the necessary condition for stability using Von Neumann is the absolute maximum of is less than or equal 1 and it is clearly satisfied then the scheme is unconditionally stable according to Von Neumann stability analysis.
5. Numerical Results
In this section we conduct some typical numerical examples to verify the accuracy, conservation laws, computational efficiency and some physical interaction phenomena described by 3-coupled nonlinear Schrödinger equations.
5.1. Single Soliton
In this test, we choose the initial condition as
The following set of parameters are used
The conserved quantities and the error for our scheme are displayed in Table 1. We have noticed that the method is conserved the conserved quantities exactly and highly accurate results are obtained. The profile of and at different times are displayed in Figure 1, Figure 2 and Figure 3 respectively.
To test the convergent rate in space and time of the proposed schemes. We define the error norm by
where and are respectively the exact and the numerical solution at the grid point . In this experiment, we take .
The convergent rate “order” is calculated by the formula.
Table 1. Errors & conserved quantities of single solitons.
Figure 1. Simulation of single soliton .
Figure 2. Simulation of single soliton .
Figure 3. Simulation of single soliton .
order (rate of convergent in space)
order (rate of convergent in time)
To calculate the convergent rate in space, we take the time step k sufficiently small and the error from temporal truncation is relatively small . From Table 2, we can easily see that the rate of convergent is 6 as we claim in this work.
To check the temporal convergent rate, we fix the spatial step h small enough so that the truncation from space is negligible such as . The results are given in Table 3 which indicate that the order is 2 as we claim in the text.
To improve the temporal accuracy of the proposed method, we use Richardson Extrapolation on the computed solution to eliminate the lower-order term in the truncation error.
Since our method applied to the scheme is in the form , we use
to eliminate the term , which makes the final solution fourth-order accurate in time dimension. Although the extrapolation requires two times as
Table 2. Rate of convergence of single solitons ( ).
Table 3. Rate of convergence of single solitons ( ).
much computation as the original scheme plus the application of the formula (47). By using the parameters
with the extrapolation formula (47), we obtain results which are displayed in Table 4.
5.2. Interaction of Two Solitons
To study the interaction of two solitons with different parameters, we choose the initial condition as a sum of two single solitons of the form
For all examples in the case of interaction, we choose the set of parameters
together with different values of for each test. We will study the dynamics of the following cases.
5.2.1. Case 1
In this test we choose the set of parameters
For this test, we have noticed that
which gives us elastic interaction. The interaction scenario is displayed in Figure 4.
Table 4. Richardson extrapolation using norm.
5.2.2. Case 2
In this test we choose the set of parameters
For this test, we have noticed that the formula
is unsatisfied which gives us inelastic interaction and it is clear in Figure 5.
For all cases, the conserved quantities given in Table 5, we have noticed that our method is conserved the conserved quantities exactly.
Figure 4. Elastic interaction of two solitons.
Figure 5. Inelastic interaction of two solitons.
Table 5. Conserved quantities of interaction of two solitons.
In this work, we have derived a highly accurate finite difference scheme for solving the 3-coupled nonlinear Schrödinger equation. The scheme is of sixth order in space and second order in time, it is unconditionally stable. A fixed point is used to solve the nonlinear block penta-diagonal system obtained. Single soliton solution and the conserved quantities are used to highlight the robustness of the method. The interaction of two solitons is discussed in detail for different parameters to produce elastic and inelastic interactions. This behavior is agreeing with    with the highest accuracy. The derived method can be used to solve similar nonlinear problems.
 Kanna, T. and Lakshmanan, M. (2003) Exact Soliton Solutions of Coupled Nonlinear Schrödinger Equation: Shape-Changing Collisions, Logic Gates and Partially Co-Herent Soliton. Physical Review, 67, 046617-046625.
 Kanna, T. and Lakshmanan, M. (2001) Exact Soliton Solutions, Shape Changing Collisions and Partially Coherent Solitons in Coupled Nonlinear Schrödinger Equations. Physical Review Letter, 86, 5043-5046.
 Lü X. and Lin, F.H. (2016) Soliton Excitation and Shape-Changing Collisions in Alpha Helical Portions with Interspine Coupling at Higher Order. Communications in Nonlinear Science and Numerical Simulation, 32, 241-261.
 Lü, X. (2014) Bright-Soliton Collisions with Shape Change by Intensity Redistribution for the Coupled Sasa-Satsuma System in the Optical Fiber Communications. Communications in Nonlinear Science and Numerical Simulation, 19, 3969-3987.
 Ismail, M.S. (2008) Numerical Solution of Coupled Nonlinear Schrödinger Equation by Galerkin Method. Mathematics and Computers in Simulation, 78, 532-547.
 Ismail, M.S., Al-Basyouni, K.S. and Aydin, A. (2015) Conservative Finite Difference Schemes for the Chiral Nonlinear Schrödinger Equation. Boundary Value Problem, 2015, Article No. 89.
 Ismaila, M.S. and Tahab, T.R. (2001) Numerical Simulation of Coupled Nonlinear Schrödinger Equation. Mathematics and Computers in Simulation, 56, 547-562.
 Ismail, M.S. and Alamri, S.Z. (2004) Highly Accurate Finite Difference Method for Coupled Nonlinear Schrödinger Equation. International Journal of Computer Mathematics, 81, 333-351.
 Ismail, M.S. and Tahab, T.R. (2007) A Linearly Implicit Conservative Scheme for the Coupled Nonlinear Schrödinger Equation. Mathematics and Computers in Simulation, 74, 302-311.
 Hu, X.L. and Zhang, L.M. (2014) Conservative Compact Difference Schemes for the Coupled Nonlinear Schrödinger System. Numerical Methods Partial Differential Equations, 30, 749-772.
 Kong, L.H., Hong, J.L., Ji, L.H. and Zhu, P.F. (2015) Compact and Efficient Conservative Schemes for Coupled Nonlinear Schrödinger Equations. Numerical Methods Partial Differential Equations, 31, 1814-1843.