Many topics in nonlinear science relating with ocean, aeronautical, mechanical, and many more may be characterized as solving nonlinear partial differential equations (NPDEs) that result from important models with mathematical and physical relevance. Solitons are one of the most frequently in the context of solutions for NPDEs, and they play a critical role in nonlinear physical phenomena exploration [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13]. In disciplines including quantum electronics, plasma physics, nonlinear optics, fluid dynamics, and many more, solitons are used to better understand the nature of nonlinear media. In addition, soliton propagation has been used to explain several important oceanic phenomena that correlate to nonlinear shallow or deep-water wave propagation [14], [15], [16], [17]. In recent years, the search for accurate soliton solutions to NPDEs has become a fascinating study topic in the field of applied sciences and engineering. To examine such NPDEs, numerous well-established techniques have been offered such as the ansatz method [18], reduction of order method [19], tanh method [20], sine-Gordon method [21], extended auxiliary function method [22], Jacobi elliptic function method [23], $F$-expansion method [24], new extended direct algebraic method [25], $(G^{\prime }/G)$-expansion method [26], homogeneous balance method [27], Lie symmetry method [28], sub-equation method [29], new Kudryashov’s method [30], simple equation method [31], Hirota bilinear method [32], and more.

One of the most general models for describing various physical nonlinear systems is the cubic nonlinear Schrödinger (CNLS) equation. This model has been used to better explain processes such as atomic physics and nonlinear optics, as well as deep sea waves, plasmas, rogue waves, and other phenomena. On the other hand, the KdV equation is one of the most prominent NPDEs. Examples include hydrodynamics, quantum field theory, plasma physics, water waves, and many more disciplines of applied sciences and engineering. It’s also utilized to explain the unidirectional propagation of long waves in nonlinear dispersive phenomena. In applied science and engineering, many types of coupled nonlinear systems have evolved as models for capturing interacting wave phenomena. Consequently, the coupled CNLS-KdV equations have been used to characterize a variety of wave phenomena in the aforementioned fields, including Langmuir waves, dust-acoustic waves, and electromagnetic waves in plasma physics, among others. They also appear in fluid mechanics phenomena such as the interactions of capillary gravity water waves, which involve short and long dispersive waves [33], [34], [35], [36], [37], [38].

In this work, the following generalized coupled CNLS-KdV equations are considered as:

where $u=u(x,t)$ represents a complex function, and $v=v(x,t)$ represents a real-valued function. The variables ${\alpha }_j$ and ${\gamma }_j,j=1,2,3$ are all real constants. In addition, $u$ denotes the short wave, and $v$ represents the long wave. In [39], [40], the existence of solutions for the Eq. (1) was highlighted. The applications of Kudryashov and sub-equation methods to the proposed system have been published in Akinyemi et al. [41]. The numerical soutions with the help of q-homotopy analysis transform method was reported in Akinyemi et al. [42]. However, our aim is to further complement our previous study on Eq. (1) by proposing a modified Sardar sub-equation (MSSE) procedure to study this coupled system. The proposed method is an innovative technique for obtaining exact solutions to differential equations. It provides a straightforward approach for dealing with nonlinear evolution equation solutions. The SSE method has been used to obtain desirable findings and assistance in the investigation of solutions to numerous problems that have arisen in applied mathematics and physics [43], [44], [45]. Other advantage of this method is that it contain solutions reported in Akinyemi et al. [41], 42] and many other new soliton type solutions. These types of solutions are essential because they can aid in the understanding of several physical processes connected to wave propagation. It is important to highlight that the results acquired in this study are completely new and has never been published before. Kumar et al. [14] used the sine-Gordon equation approach to investigate the generalized Schrödinger–Boussinesq equations, which is another coupled model that describes the interaction between complex short and real long waves envelope.

The following is how the remainder of this paper is structured: Section 2 presents the description and application of recommended techniques. In Section 3, we will study several types of solitions type solutions of Eq. (1) with their valid constraint conditions, as well as discuss their dynamic characteristics using graphs, utilizing the MSSE method. Finally, in Section 5, the conclusion is provided.

This method is a powerful and robust approach that may be used to generate numerous types of soliton solutions, including dark, bright, W-shaped, mixed dark-bright, singular, mixed singular solitons, periodic, and other solutions, for both the classical and fractional order NPDEs as compared to other methods in the literature. The fact that this approach overcomes the complications of the solitary wave ansatz method [18], [46] is noteworthy. As a result, this section contains a full description of the MSSE approach. To use this approach, we assume the following general NPDE as:

where $G$ denotes a function of $u$ and its derivatives. The complex wave transformation [41] given as

is proposed to reduce Eq. (2) into nonlinear ordinary differential equations (NODE):

Here, $U=U\left(\zeta \right)$ denotes the amplitudes while ${\omega }_j$ and ${\eta }_j,j=1,2$ are arbitrary constants. The solutions of Eq. (4) can be classify as:

with the constants ${\delta }_r,(r=0,1,\dots ,M),$ ${\omega }_j,$ and ${\eta }_j,(j=1,2)$are to be computed later. Balancing the highest-order derivative and the nonlinear terms of Eq. (4), determines the integer $M.$ The function $\mathcal{Q}\left(\zeta \right)$ in Eq. (5) satisfies

where ${\vartheta }_0,{\vartheta }_1,$ and ${\vartheta }_2$are constants. In addition, the family of solutions for Eq. (6) with constant $\varrho$ are listed as below:

1. If ${\vartheta }_0=0,{\vartheta }_1>0,$ and ${\vartheta }_2\ne 0,$ then

2. For constants $A_1$ and $A_2.$ Let ${\vartheta }_2=\pm 4A_1{A}_2,{\vartheta }_1>0,$ and ${\vartheta }_0=0,$ we have

3. If ${\vartheta }_0=\frac{{\vartheta }_1^2}{4{\vartheta }_2},{\vartheta }_1<0,{\vartheta }_2>0,$ with constants $B_1,$ and $B_2,$ then

4. If ${\vartheta }_0=0,{\vartheta }_1<0,$ and ${\vartheta }_2\ne 0,$ then

5. If ${\vartheta }_0=\frac{{\vartheta }_1^2}{4{\vartheta }_2},{\vartheta }_1>0,{\vartheta }_2>0,$ and $B_1^2-B_2^2>0,$ then

6. If ${\vartheta }_0=0$ and ${\vartheta }_1>0,$ then

7. If ${\vartheta }_0={\vartheta }_1=0$ and ${\vartheta }_2>0,$ then

8. If ${\vartheta }_0={\vartheta }_1=0$ and ${\vartheta }_2<0,$ then

Inserting Eqs. (5) and (6) into Eq. (4), then gathering and equating all the coefficients of ${\mathcal{Q}}^r\left(\zeta \right)$ to zero yields an algebraic system of equations in unknown variables ${\delta }^r,{\omega }_j,$ and ${\eta }_j.$ After that, we solve the relevant set of equations, place the variables in Eq. (5), and use Eqs. (7)–(27) to obtain the precise solutions of Eq. (2).

As presented in Eq. (1), the CNLS-KdV equations as a generalized coupled system:

Because $u$ is a complex function and $v$ is a real-valued function, we suggest the transformation as follows:

where ${\omega }_j$ and ${\eta }_j,j=1,2$ are arbitrary constants. Using the transformation defined in Eq. (29), we obtain the real and imaginary parts of Eq. (28) as follows:

and

Solving Eq. (31) yields

Inserting Eq. (33) into Eq. (32), then integrating once while disregarding the integration constant results

From Eq. (5), we can express the solutions of Eqs. (30) and (34) as:

Here, $f_r$ and $g_r$ are constant to be evaluated later. Balancing $U^{″}$ and $U^3$ in Eq. (30) results $M_1=1,$ while $V^{″}$ and $V^2$ in Eq. (34) yields $M_2=2.$

The solutions take the form

Putting Eqs. (6) and (36) into Eqs. (30) and (34), gathering all the coefficients of ${\mathcal{Q}}^r\left(\zeta \right)$ and setting them to zero, we get

The solutions to the aforementioned algebraic equations using Mathematica, along with Eq. (33) yield the following results:

By incorporating the parameters in Case 1 into Eq. (36), in addition to Eqs. (7)–(27), we have the following solutions:

The bright and singular solitons can be consider as:

where $\zeta ={\omega }_{1}x+{\eta }_{1}t,{\vartheta }_0=0,{\vartheta }_1>0,$ and ${\vartheta }_2\ne 0.$

The dark and singular solitons can be consider as:

where ${\vartheta }_0=\frac{{\vartheta }_1^2}{4{\vartheta }_2},{\vartheta }_1<0,{\vartheta }_2>0,$ and $H_1=-\frac{6{\gamma }_2{\gamma }_1(4{\gamma }_2{\vartheta }_1{\alpha }_2{\omega }_1^2-{\gamma }_3{\alpha }_3)+\sqrt{{\gamma }_1^2((24{\gamma }_2^2{\vartheta }_1{\alpha }_2{\omega }_1^2+{\gamma }_3({\gamma }_1{\alpha }_1-6{\gamma }_2{\alpha }_3)){}^2-432{\gamma }_2^4{\vartheta }_1^2{\alpha }_2^2{\omega }_1^4)}+{\gamma }_3{\gamma }_1^2{\alpha }_1}{6{\gamma }_1^2{\gamma }_2{\alpha }_2}.$

The combo bright-singular solitons are expressed as:

where ${\vartheta }_0=0,{\vartheta }_1>0,{\vartheta }_2=\pm 4A_1{A}_2,$ and $H_2=-\frac{6{\gamma }_2{\gamma }_1(4{\gamma }_2{\vartheta }_1{\alpha }_2{\omega }_1^2-{\gamma }_3{\alpha }_3)+\sqrt{{\gamma }_1^2(24{\gamma }_2^2{\vartheta }_1{\alpha }_2{\omega }_1^2+{\gamma }_3({\gamma }_1{\alpha }_1-6{\gamma }_2{\alpha }_3)){}^2}+{\gamma }_3{\gamma }_1^2{\alpha }_1}{6{\gamma }_1^2{\gamma }_2{\alpha }_2}.$

The combo solitons are as follows:

where ${\vartheta }_0=\frac{{\vartheta }_1^2}{4{\vartheta }_2},{\vartheta }_1<0,{\vartheta }_2>0,$ and $H_1=-\frac{6{\gamma }_2{\gamma }_1(4{\gamma }_2{\vartheta }_1{\alpha }_2{\omega }_1^2-{\gamma }_3{\alpha }_3)+\sqrt{{\gamma }_1^2((24{\gamma }_2^2{\vartheta }_1{\alpha }_2{\omega }_1^2+{\gamma }_3({\gamma }_1{\alpha }_1-6{\gamma }_2{\alpha }_3)){}^2-432{\gamma }_2^4{\vartheta }_1^2{\alpha }_2^2{\omega }_1^4)}+{\gamma }_3{\gamma }_1^2{\alpha }_1}{6{\gamma }_1^2{\gamma }_2{\alpha }_2}.$

The trigonometric function solutions are acquired as follows:

for ${\vartheta }_0=0,{\vartheta }_1<0,$ and ${\vartheta }_2\ne 0,$ and $H_2=-\frac{6{\gamma }_2{\gamma }_1(4{\gamma }_2{\vartheta }_1{\alpha }_2{\omega }_1^2-{\gamma }_3{\alpha }_3)+\sqrt{{\gamma }_1^2(24{\gamma }_2^2{\vartheta }_1{\alpha }_2{\omega }_1^2+{\gamma }_3({\gamma }_1{\alpha }_1-6{\gamma }_2{\alpha }_3)){}^2}+{\gamma }_3{\gamma }_1^2{\alpha }_1}{6{\gamma }_1^2{\gamma }_2{\alpha }_2}.$

for ${\vartheta }_0=\frac{{\vartheta }_1^2}{4{\vartheta }_2},{\vartheta }_1>0,{\vartheta }_2>0,$ and $H_1=-\frac{6{\gamma }_2{\gamma }_1(4{\gamma }_2{\vartheta }_1{\alpha }_2{\omega }_1^2-{\gamma }_3{\alpha }_3)+\sqrt{{\gamma }_1^2((24{\gamma }_2^2{\vartheta }_1{\alpha }_2{\omega }_1^2+{\gamma }_3({\gamma }_1{\alpha }_1-6{\gamma }_2{\alpha }_3)){}^2-432{\gamma }_2^4{\vartheta }_1^2{\alpha }_2^2{\omega }_1^4)}+{\gamma }_3{\gamma }_1^2{\alpha }_1}{6{\gamma }_1^2{\gamma }_2{\alpha }_2}.$

The mixed trigonometric function solutions are as follows: