EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6490 ISSN 1307-5543 – ejpam.com Published by New York Business Global Innovative Solitary Wave Solutions for the (3+1)-Dimensional Boussinesq Kadomtsev-Petviashvili-Type Equation Derived via the Improved Modified Extended Tanh-Function Method Wael W. Mohammed1, Abeer S. Khalifa2, Hijyah Alshammary1, Hamdy M. Ahmed3, Mohamed S. Algolam1, Reda Elbarougy4, Karim K. Ahmed5,∗ 1 Department of Mathematics, College of Science, University of Ha’il, Ha’il 2440, Saudi Arabia 2 Department of Mathematics, Faculty of Basic Sciences, The German University in Cairo (GUC), Cairo, Egypt 3 Department of Physics and Engineering Mathematics, Higher Institute of Engineering, El Shorouk Academy, Cairo, Egypt 4 Department of Artificial Intelligence and Data Science, College of Computer Science and Engineering, University of Ha’il, Saudi Arabia 5 Department of Mathematics, Faculty of Engineering, German International University (GIU), New Administrative Capital, Cairo, Egypt Abstract. The goal of this research is to create and analyze a novel (3+1) dimensional model that incorporates two different equations: a three-dimensional Kadomtsev-Petviashvili equation and a three-dimensional Boussinesq-KP-type equation. One of the unexpected outcomes of the idea of mixing integrable equations is a resonance of solitons. This paper presents a wide range of possible analytical solutions for the pKP–BKP equation in (3+1)-dimensions, including dark, bright, singular solitons, and other exact solutions like singular periodic, Jacobi elliptic function, rational, and exponential type. The (3+1)-dimensional B-KP-type model is subjected to the im- proved modified extended tanh-function approach in order to obtain novel traveling wave solutions. The employed equation plays a crucial role in describing and interpreting a broad range of non- linear phenomena seen in fluid mechanics and other nonlinear engineering and physics issues due to the strong correlation and wide range of applications of the Boussinesq-type and KP equations. The approach can help to find other kinds of solutions to the chosen equation that have not been found and published in the literature before. These solutions can aid in the comprehension of wave propagation in water wave dynamics. To further facilitate learning, they are replicated through the use of contour graphics, 2D, and 3D symbolic calculations. Moreover, linear stability analysis is discussed for the obtained solutions. 2020 Mathematics Subject Classifications: 35C05, 35C07, 35C08 Key Words and Phrases: Soliton solutions, Boussinesq-type equation, Kadomtsev-Petviashvili equation, NPDEs, linear stability analysis ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6490 Email addresses: karim.kamal.502@gmail.com (K. K. Ahmed) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 2 of 23 1. Introduction A variety of nonlinear evolution equations, such as Boussinesq, Kadomtsev–Petviashvili, and nonlinear Schrödinger-type models, are examined in some latest research, with an em- phasis on building precise soliton and wave solutions by analytical techniques [1–5]. Mod- eling nonlinear wave propagation, optical solitons, and interaction dynamics in a variety of physical media, including fluids, birefringent fibers, magneto-optic waveguides, and twin- core couplers, is made easier by the results [6–10]. Many recent researchers have committed their efforts to examining the integrability of various nonlinear evolution equations using a variety of methods and techniques because they recognize the importance and effectiveness of studying integrability prior to analyzing distinct nonlinear evolution equations and un- derstanding their unique behavior. For example, in [11], the Painlevé property for partial differential equations (PDEs) was explained in detail, and its significance in determining the integrability was illustrated. The Painlevé property for the PDEs was explained in detail, and its significance in determining the integrability was illustrated. The Painlevé test was utilized by the authors to verify the integrability of numerous universal evolu- tion equations, including Burgers, Sine-Gordon, Boussinesq (B-type), Korteweg-de Vries (KdV)-type, and the Kadomtsev-Petviashvili (KP) equations. In many scientific domains, these equations may be used to mimic a wide range of nonlinear processes. Several studies in the literature have been effective in providing an explanation for a wide range of non- linear and enigmatic phenomena that appear in various physical systems due to numerous different evolution equations that have been studied for integrability. For instance, in [12– 14], the properties of either modulated or unmodulated solitary waves (SWs) and cnoidal waves (CWs) that exist in fluid mechanics, seas, and oceans, as well as in narrow channels, have been effectively investigated using the KdV-type equation and other equations (such as the nonlinear Schrödinger equation (NLSE) and Boussinesq-type equations). These formulas have also been used to investigate nonlinear optical communications, solitons, and shocks in plasma physics. It is established that in the system under study, a balance between nonlinearity and dispersion can produce solitons. The SWs become shock waves when the dispersion and nonlinearity of the system are unbalanced. Assume that the system has a dissipation force, dispersion, and nonlinearity, and that the dissipation force predominates over the dispersion. The solitons will then change into shock waves in such scenarios. Recent studies have applied a host of analysis techniques to obtain exact solutions of higher-dimensional and nonlinear fractional evolution equations that are significant in wave dynamics and nonlinear optics. For instance, Almatrafi [15, 16] applied the improved modified extended tanh-function method and other analytical techniques to obtain solitary wave solutions of fractional models with emphasis on lower-dimensional equations with- out taking into account the additional complexity arising from higher spatial dimensions. Murad et al. [17] examined the Sasa–Satsuma equation of higher-order time-fractional in optical fibers, and Younas et al. [18] considered interaction phenomena in the (2+1)- dimensional KdV–Sawada–Kotera–Ramani equation. Younas et al. [19, 20] also discussed soliton dynamics in optical systems and magneto-electro-elastic media. However, these W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 3 of 23 studies were concentrated on specific physical models with a predetermined dimension- ality or special nonlinear structures. On the other hand, the present work is concerned with the dynamics of a novel (3+1)-dimensional nonlinear wave equation with conformable fractional derivatives, in which a more generalized model frame is presented that can por- tray complex nonlinear features in higher-dimensional fluid and optical systems. This extension fills the existing gap in the literature by providing a more comprehensive class of exact solutions, including bright, dark, singular, and periodic solitons, and thus devel- oping additional insight into the spatiotemporal dynamics of nonlinear waves in realistic multidimensional settings. The Boussinesq equation (BE), developed by Boussinesq in 1871, is one of the most important evolution equations used in plasma physics and fluid mechanics. The propaga- tion of small-amplitude shallow water waves traveling at a constant speed in a continuous depth water channel is described by this equation and certain extensions from it [21]. Numerous studies have tackled this equation, scrutinizing and resolving it using diverse methodologies. For example, Wazwaz [22] used a family of tanh methods to evaluate it and get only a periodic solution as well as a single soliton solution. To obtain multiple soliton solutions for this model, the author also used Hirota’s direct approach in conjunc- tion with its simplified form. In [23], the author obtained solutions in the form of single solitons and periodic waves by analyzing the behavior of this equation using the modified Adomian decomposition technique. Furthermore, a variety of nonlinear phenomena seen in a wide range of engineering and physical systems have been accurately modeled by the KdV equations. The KP equation (KPE) belongs to the KdV family of evolutionary equations, partic- ularly when two-dimensional propagation and perturbation are taken into account as in [24]. In this study, the two-dimensional KdV equation has an unlimited number of con- served quantities and is fully integrable, which relates to the simulation of many nonlinear structures in a wide range of real-world circumstances, including a harmonic lattice, a two-layer liquid with gradually varying depth, and plasma physics. Inspired by several applications of both BE and KPE, the goal of our present work is to combine the two equations to produce a new model that may satisfactorily describe a variety of occurrences for which neither equation could account. As in [25], the work in- troduced a brand-new, three-dimensional Boussinesq-KP-type (B-KP-type) equation with four linear terms that were balanced with a higher order linear term of fourth order that represented the dispersion effect. This model was analyzed using the tanh-method family; however, periodic wave and single soliton solutions were only obtained. In this study, the improved modified extended (IME) tanh-function method is sug- gested as one simple method to analyze the novel constructed evolution equations in [25], which is named as (3+1)-dimensional B-KP-type equation. Several traveling wave and soliton solutions are derived and discussed. A wide range of complex and nonlinear phe- nomena that develop and spread in different nonlinear engineering and physical systems, particularly the nonlinear phenomena found in fluid dynamics, were clearly explained through some 2D, 3D, and contour illustrations that show propagation properties and behaviors of some obtained traveling wave solutions. In addition, to discuss the stability W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 4 of 23 of the obtained solutions, we present linear stability analysis in some cases, like neutrally stable and instability. The (3+1)-dimensional B-KP-type equation is given as follows: Cϕtt + ϕxxxx + αϕxx + βϕxy + γϕxz + ϑϕxt + ϱ(ϕ2)xx + µϕyy = 0, (1) where ϕ ≡ ϕ(x, y, z, t) represents a real-valued function that must be sufficiently frequently differentiable, and subscripts denote partial derivatives. It refers to the height of a fluid’s free surface; x, y, z refer to the variables of the space dimensions, while t denotes the time variation. C, α, β, γ, ϑ, ϱ and µ are arbitrary real values parameters that will be calculated later during this study. ϕ appears with its partial derivatives to describe the fluid wave profiles with significant amplitudes that are maintained for a brief duration of time, which signifies localization in both the space and the time domains. The (3+1)- dimensional pKP-BKP equation, which depicts the propagation of nonlinear waves in three spatial dimensions and one temporal dimension, is particularly significant for studying complex wave dynamics in systems like fluids, plasmas, and optical media. It extends the standard KP equation by including wave interactions that occur both longitudinally and transversely, and it shows how perturbations in one direction can propagate over many spatial dimensions. Understanding this equation is necessary for understanding phenomena such as shallow water waves, ion-acoustic waves in plasmas, and even optical pulses in nonlinear media, where the stability and growth of the wave crucially depend on the balance between nonlinearity and dispersion. The main objective of this approach is to maximize the use of Eq. (1) by including its parameters to generate more types of analytical solutions when utilized in an appropriate, straightforward, easy manner, and with the assistance of symbolic computations. Numer- ous studies have employed the IME tanh-function approach to examine a broad spectrum of solitons and other precise wave solutions [13]. This study uses the IME tanh-function approach for Eq. (1) to compute various solitons and other exact wave solutions. This particular combination is novel research that has never been done before, and it uses the IME tanh-function approach to provide accurate analytical solutions to nonlinear partial differential equations (NPDEs) that are insightful, efficient, and versatile. Numerous so- lutions are obtained using this approach, such as the Jacobi epsilon function, exponential, singular periodic, rational, dark soliton, and bright soliton solutions. Furthermore, the recovered solutions confirm the effectiveness and strength of the used approach. The structure of this article is as follows. An outline of the suggested model and its theoretical underpinnings is given in Section 1. The main components of the IME tanh- function algorithm are presented in Section 2. In Section 3, the symbolic computations are completed and the results are summarized using the Wolfram Mathematica program. Section 4 presents multiple dynamic wave patterns of different soliton solutions graphically using 3-D, contour, and 2-D simulations and their discussion. Section 5 presents a novel comparison with the most common methods in the literature. An implementation of the linear stability analysis on the governing model is discussed in Section 6. The study conclusions are reported in Section 7. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 5 of 23 2. Improved modified extended tanh-function method The general procedures of the IME tanh-function method [26, 27] will be shown in this section, in addition to its motivation and advantages in solving NPDE. 2.1. Quick view of the method As we start looking at the NPDE that’s shown below [28, 29]: J (ϕ, ϕx, ϕy, ϕz, ϕt, ϕxx, ϕxy, ϕxz, . . .) = 0, (2) where J is a polynomial that is composed of ϕ(x, y, z, t) and some of ϕ partial derivatives with respect to both time (t) and the space dimensions of the used dynamic system (x, y, z). Procedure-(I): Considering the following wave transformation ϕ(x, y, z, t) = R(η); η = x+ ky + ℵz − Ωt, (3) in which k , ℵ and Ω are real valued constants. R acts as the function of the obtained solution. Considering Eq. (3) and Eq. (2) together, with performing rearranging, the following nonlinear ordinary differential equation (NLODE) form is derived: Z(R, R′, R′′, R′′′, . . .) = 0. (4) Procedure-(II): In order to achieve the solution of Eq. (4) according to the applied method, the solution is proposed in the following truncated series: R(η) = N∑ i=0 diW i(η) + N∑ i=1 fiW−i(η), (5) where d0, d1 ..., dN and f1, ...,fN are real constants to be calculated, under condition that dN and fN should not be zero, simultaneously. Procedure-(III): By applying the homogeneous balance principle between the nonlin- earity and the dispersion of Eq. (4) to determine the value of the balancing constant N. In addition to considering the following constraint for W(η): W ′(η) = ϵ √ τ0 + τ1W(η) + τ2W2(η) + τ3W3(η) + τ4W4(η), (6) where ϵ = ±1 and τm (0 ≤ m ≤ 4) are real-valued constants. Several types of fundamental solutions are obtained from Eq. (6) using the many potential values of τ0, τ1, τ2, τ3 and τ4 which are shown in Appendix A. Procedure-(IV): An equation in powers of W(η) is obtained by inserting Eqs. (5) and (6) into Eq. (4). The system of algebraic non-linear equations that results from applying algebraic polynomial operations on the coefficients of W(η) and equating them to zero may be solved with a variety of software applications, including the Wolfram Mathematica. This will finally lead to the production of several exact solutions for Eq. (2). The flowchart in figure (1) shows a brief description of the used scheme. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 6 of 23 Figure 1: Flowchart of the IME tanh-function method. 2.2. Motivation and advantages of the method The IME tanh-function method is highly advantageous for obtaining analytical so- lutions and gaining a deeper understanding of soliton dynamics when compared to the other most current approaches. However, it might be challenging to regulate because of its complexity and sensitivity to the initial conditions. Depending on the specifics of the model being used, such as the type of NPDEs, determining the boundary constraints, and the desired ratio of computational practicality to analytical proficiency, each technique is appropriate. Alternative techniques are more complex, but they also allow for greater flexibility and improvement. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 7 of 23 3. Extraction of solitons and some other solutions By applying Eq. (3) for Eq. (1), then Eq. (1) appears as a non-linear ordinary differential equation (NODE), having the following form: 2ϱ ( R′)2 + ( α+ k(β + kµ) + 2ϱR− ϑΩ+ CΩ2 + γℵ ) R′′ +R(4) = 0. (7) Thus, the exact solutions for Eq. (7) can be generated in the following manner by using the principle of balance that was mentioned in Section 2 to Eq. (7), the balance is performed between the terms include R(4) and RR′′, determining the balance constant N = 2. Hence, Eq. (5) becomes R(η) = d0 + d1W(η) + d2W2(η) + f1 W(η) + f2 W2(η) . (8) Using Eq. (8) and the constraint in Eq. (6), Eq. (7) yields a polynomial in W(η). An algebraic non-linear equations system is created when all terms with the same power are combined and set to equal zero. The Wolfram Mathematica software tool is used to solve these equations and generate the potential solutions. It is specified that d2 and f2 cannot both be zero simultaneously. Family (1): If τ0 = τ1 = τ3 = 0, the below set of solutions is resulted: d1 = f1 = f2 = 0, d0 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ , d2 = −6τ4 ϱ . The following exact solutions are obtained for Eq. (1), as per the above set of solutions: (1.1) If τ2 > 0, τ4 < 0 and ϱ ̸= 0, the following bright soliton solution is obatined: ϕ1.1 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6τ2 ϱ sech2 [(x+ ky + ℵz − Ωt) √ τ2] . (9) (1.2) If τ2 < 0, τ4 > 0 and ϱ ̸= 0, the following singular periodic solution is obtained: ϕ1.2 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6τ2 ϱ sec2 [ (x+ ky + ℵz − Ωt) √ −τ2 ] . (10) (1.3) If τ2 = 0 and τ4 > 0, the following rational solution; (ϱ)(x+ ky + ℵz − Ωt) ̸= 0: ϕ1.3 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) 2ϱ − 6 ϱ(x+ ky + ℵz − Ωt)2 . (11) Family (2): If τ1 = τ3 = 0, the sets of solutions listed below are: (2.1) d1 = f1 = f2 = 0, d0 = −α+γℵ−Ω(ϑ−CΩ)+k(β+kµ)+4τ2 2ϱ , d2 = −6τ4 ϱ . W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 8 of 23 (2.2) d1 = d2 = f1 = 0, d0 = −α+γℵ−Ω(ϑ−CΩ)+k(β+kµ)+4τ2 2ϱ , f2 = −6τ0 ϱ . (2.3) d1 = f1 = 0, d0 = −α+γℵ−Ω(ϑ−CΩ)+k(β+kµ)+4τ2 2ϱ , d2 = −6τ4 ϱ , f2 = −6τ0 ϱ . By considering the solutions’ set (2.1), the corresponding exact solutions of Eq. (1) are produced as: (2.1,1) If τ0 = τ22 4τ4 , τ2 < 0, τ4 > 0 and ϱ ̸= 0, the following dark soliton solution is obtained: ϕ2.1,1 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 3τ2 ϱ tanh2 [ (x+ ky + ℵz − Ωt) √ −τ2 2 ] . (12) (2.1,2) If τ0 = τ22 4τ4 , τ2 > 0, τ4 > 0 and ϱ ̸= 0, the following singular periodic solution is formed: ϕ2.1,2 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ −3τ2 ϱ tan2 [ (x+ ky + ℵz − Ωt) √ τ2 2 ] . (13) (2.1,3) If τ0 = m2(1−m2)τ22 (2m2−1)2τ4 , τ2 > 0, τ4 < 0, 1√ 2 < m ≤ 1 and ϱ ̸= 0, then the below Jacobi elliptic function solution (JEFS) is obtained: ϕ2.1,3 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6m2τ2 cn2 [ (x+ ky + ℵz − Ωt) √ τ2 2m2−1 ] ϱ(2m2 − 1) . (14) Set m = 1 in Eq. (14) as a special case, a bright soliton solution is produced as follows: ϕ2.1,4 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6τ2 ϱ sech2 [(x+ ky + ℵz − Ωt) √ τ2] . (15) (2.1,4) If τ0 = (1−m2)τ22 (2−m2)2τ4 , τ2 > 0, τ4 < 0, 0 < m ≤ 1 and ϱ ̸= 0, the below JEFS is raised: ϕ2.1,5 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6m2 dn2 [ (x+ ky + ℵz − Ωt) √ τ2 2−m2 ] ϱ(2−m2) . (16) Set m = 1 in Eq. (16) as a special case, its solution resulted as a bright soliton in the following form: ϕ2.1,6 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6 ϱ sech2 [(x+ ky + ℵz − Ωt) √ τ2] . (17) W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 9 of 23 (2.1,5) If τ0 = m2τ22 (m2+1)2τ4 , τ2 < 0, τ4 > 0, 0 < m ≤ 1 and ϱ ̸= 0, then a JEFS is resulted as below: ϕ2.1,7 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6m2 sn2 [ (x+ ky + ℵz − Ωt) √ − τ2 m2+1 ] ϱ(m2 + 1) . (18) Special case, when setting m = 1 in Eq. (18), the following dark soliton solution can be obtained: ϕ2.1,8 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 3τ2 ϱ tanh2 [ (x+ ky + ℵz − Ωt) √ −τ2 2 ] . (19) The mentioned set of solutions (2.2) indicates that Eq. (1) has exact solutions, which can be phrased as: (2.2,1) If τ0 = τ22 4τ4 , τ2 < 0, τ4 > 0 and ϱ ̸= 0, the below singular soliton form is appeared in the obtained solution: ϕ2.2,1 = −α+ β + kµ+ γℵ − Ω(ϑ− CΩ) + 4τ2 2ϱ + 3τ2 ϱ coth2 [ (x+ ky + ℵz − Ωt) √ −τ2 2 ] . (20) (2.2,2) If τ0 = τ22 4τ4 , τ2 > 0, τ4 > 0 and ϱ ̸= 0, following singular periodic solution is obtained: ϕ2.2,2 = −α+ β + kµ+ γℵ − Ω(ϑ− CΩ) + 4τ2 2ϱ −3τ2 ϱ cot2 [ (x+ ky + ℵz − Ωt) √ τ2 2 ] . (21) (2.2,3) If τ0 = m2(1−m2)τ22 (2m2−1)2τ4 , τ2 > 0, τ4 < 0, 1√ 2 < m < 1 and ϱ ̸= 0, the below JEFS is obtained: ϕ2.2,3 = −α+ β + kµ+ γℵ − Ω(ϑ− CΩ) 2ϱ −2τ2 ϱ 1 + 3 ( m2 − 1 ) (2m2 − 1) cn2 [ (x+ ky + ℵz − Ωt) √ τ2 2m2−1 ]  . (22) (2.2,4) If τ0 = (1−m2)τ22 (2−m2)2τ4 , τ2 > 0, τ4 < 0, 0 < m < 1 and ϱ ̸= 0, the below JEFS is raised: ϕ2.2,4 = −α+ β + kµ+ γℵ − Ω(ϑ− CΩ) + 4τ2 2ϱ + 6 ( m2 − 1 ) τ22 m2 (m2 − 2) ϱ dn2 [ (x+ ky + ℵz − Ωt) √ τ2 2−m2 ] . (23) W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 10 of 23 (2.2,5) If τ0 = m2τ22 (m2+1)2τ4 , τ2 < 0, τ4 > 0, 0 ≤ m ≤ 1 and ϱ ̸= 0, the obtained solution is resulted as the below JEFS: ϕ2.2,5 = −α+ β + kµ+ γℵ − Ω(ϑ− CΩ) 2ϱ −2τ2 ϱ 1− 3 (1 +m2) sn2 [ (x+ ky + ℵz − Ωt) √ − τ2 1+m2 ]  . (24) Set m = 0 in Eq. (24) as a special case, then the solution resulted as the below singular periodic solution: ϕ2.2,6 = −α+ β + kµ+ γℵ − Ω(ϑ− CΩ) + 4τ2 2ϱ + 6τ2 ϱ csc2 [ (x+ ky + ℵz − Ωt) √ −τ2 ] . (25) Set m = 1 in Eq. (24) as a special case, then the solution resulted as the following singular soliton solution: ϕ2.2,7 = −α+ β + kµ+ γℵ − Ω(ϑ− CΩ) + 4τ2 2ϱ + 3τ2 ϱ coth2 [ (x+ ky + ℵz − Ωt) √ −τ2 2 ] . (26) Considering the set of solutions (2.3), then Eq. (1) has exact solutions, which can be phrased as: (2.3,1) If τ0 = τ22 4τ4 , τ2 < 0, τ4 > 0, and ϱ ̸= 0, the solution is obtained in a singular soliton form given as follows: ϕ2.3,1 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 16τ2 2ϱ + 12τ2 ϱ coth2 [ (x+ ky + ℵz − Ωt) √ −2τ2 ] . (27) (2.3,2) If τ0 = τ22 4τ4 , τ2 > 0, τ4 > 0 and ϱ ̸= 0, the following singular periodic solution is obtained: ϕ2.3,2 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ)− 8τ2 2ϱ −12τ2 ϱ csc2 [ (x+ ky + ℵz − Ωt) √ 2τ2 ] . (28) (2.3,3) If τ0 = m2(1−m2)τ22 (2m2−1)2τ4 , τ2 > 0, τ4 < 0, 1√ 2 < m ≤ 1 and ϱ ̸= 0, a JEFS is obtained as follows: ϕ2.3,3 = − α + γℵ − Ω(ϑ − CΩ) + k(β + kµ) + 4τ2 2ϱ − 6τ2 ϱ  (m2 − 1)( 2m2 − 1 ) cn2 [ (x + ky + ℵz − Ωt) √ τ2 2m2−1 ] − m2cn2 [ (x + ky + ℵz − Ωt) √ τ2 2m2−1 ] 2m2 − 1  . (29) Set m = 1 in Eq. (29) as a special case, its solution is obtained as the below bright soliton solution: ϕ2.3,4 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6τ2 ϱ sech2 [(x+ ky + ℵz − Ωt) √ τ2] . (30) W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 11 of 23 (2.3,4) If τ0 = (1−m2)τ22 (2−m2)2τ4 , τ2 > 0, τ4 < 0, 0 < m ≤ 1 and ϱ ̸= 0, a JEFS is obtained as below: ϕ2.3,5 = − α + γℵ − Ω(ϑ − CΩ) + k(β + kµ) + 4τ2 2ϱ + 6 ϱ m2dn2 [ (x + ky + ℵz − Ωt) √ τ2 2−m2 ] 2 − m2 − ( m2 − 1 ) τ2 2 m2 ( 2 − m2 ) dn2 [ (x + ky + ℵz − Ωt) √ τ2 2−m2 ]  . (31) Set m = 1 in Eq. (31) as a special case, the solution is obtained as a bright soliton type as below: ϕ2.3,6 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6 ϱ sech2 [(x+ ky + ℵz − Ωt) √ τ2] . (32) (2.3,5) If τ0 = m2τ22 (m2+1)2τ4 , τ2 < 0, τ4 > 0, ϱ ̸= 0 and 0 ≤ m ≤ 1, the following JEFS is resulted: ϕ2.3,7 = − α + γℵ − Ω(ϑ − CΩ) + k(β + kµ) + 4τ2 2ϱ + 6τ2 ϱ  1( m2 + 1 ) sn2 [ (x + ky + ℵz − Ωt) √ − τ2 m2+1 ] + m2sn2 ( x + ky + ℵz − Ωt) √ − τ2 m2+1 ] m2 + 1  . (33) Special case, when setting m = 0 in Eq. (33), the following singular periodic solution can be obtained: ϕ2.3,8 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 4τ2 2ϱ + 6τ2 ϱ csc2 [ (x+ ky + ℵz − Ωt) √ −τ2 ] . (34) Special case, when setting m = 1 in Eq. (33), the following singular soliton solution can be obtained: ϕ2.3,9 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ) + 16τ2 2ϱ + 12τ2 ϱ coth2 [ (x+ ky + ℵz − Ωt) √ −2τ2 ] . (35) Family (3): If τ2 = τ4 = 0, the resulted set of solutions is mentioned below: d1 = d2 = 0, d0 = −α+ γℵ − Ω(ϑ− CΩ) + k(β + kµ)− 3 3 √ τ0τ23 2ϱ , f1 = 6 3 √ τ20 τ3 ϱ , f2 = −6τ0 ϱ , τ1 = −2 3 √ τ20 τ3. The exact solution to Eq. (1) that arises from this set of solutions is as follows: If τ0 ̸= 0, τ1 ̸= 0, τ3 > 0 and ϱ ̸= 0, a Weierstrass elliptic doubly periodic function is produced in the form of solution: ϕ3.1 = − α − Ω(ϑ − CΩ) + k(β + kµ) − 3 3 √ τ0τ 2 3 2ϱ + 6 ( τ0 − 3 √ τ2 0 τ3 ℘ [ (x + ky + ℵz − Ωt) √ τ3 4 ; ( 8 3 √ τ2 0 τ3 τ3 ,− 4τ0 τ3 ) ]) ϱ ℘2 [ (x + ky + ℵz − Ωt) √ τ3 4 ; ( 8 3 √ τ2 0 τ3 τ3 ,− 4τ0 τ3 ) ] . (36) Family (4): If τ3 = τ4 = 0, the raised sets of solutions are generated as: (4.1) d1 = d2 = 0, d0 = −α+γℵ−Ω(CΩ+ϑ)+k(β+kµ)+τ2 2ϱ , f1 = 6 √ τ0τ2 ϱ , f2 = − 6τ0 ϱ , τ0 = τ2 1 4τ2 . W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 12 of 23 (4.2) d1 = f1 = d2 = τ1 = 0, d0 = −α+γℵ−Ω(CΩ+ϑ)+k(β+kµ)+4τ2 2ϱ , f2 = − 6τ0 ϱ . Using the revealed set of solutions (4.1), Eq. (1) gives its exact solution that can be expressed as: (4.1) If τ0 > 0, τ1 > 0, τ2 > 0 and ϱ ̸= 0, the following exponential solution is produced such that τ1 − 2τ2e √ τ2(x+ky+ℵz−Ωt) ̸= 0: ϕ4.1 = −α+ γℵ − Ω(CΩ+ ϑ) + k(β + kµ) + τ2 2ϱ − 12τ1τ2 ϱ [ τ1 − τ2e (x+ky+ℵz−Ωt) √ τ2( τ1 − 2τ2e(x+ky+ℵz−Ωt) √ τ2 )2 ] . (37) By inserting the above set of solutions (4.2) to Eq. (1), the following solutions are obtained: (4.2,1) If τ0 > 0, τ2 < 0 and ϱ ̸= 0, a singular periodic solution is produced in the following form: ϕ4.2,1 = −α+ γℵ − Ω(CΩ+ ϑ) + k(β + kµ) + 4τ2 2ϱ − 6τ2 ϱ csc2 [ (x+ ky + ℵz − Ωt) √ −τ2 ] . (38) (4.2,2) If τ0 > 0, τ2 > 0 and ϱ ̸= 0, a singular soliton solution is produced in the following form: ϕ4.2,2 = −α+ γℵ − Ω(CΩ+ ϑ) + k(β + kµ) + 4τ2 2ϱ − 6τ2 ϱ csch2 [(x+ ky + ℵz − Ωt) √ τ2] . (39) Family (5): If τ0 = τ1 = d0 = 0 and τ4 > 0, the sets of the resulted solutions are listed below (5.1) d1 = f1 = f2 = τ3 = 0, d2 = − 6τ4 ϱ , ϑ = α+γℵ+k(β+kµ)+CΩ2+4τ2 Ω . (5.2) f1 = f2 = 0, d1 = 6 √ τ2τ4 ϱ , d2 = − 6τ4 ϱ , τ2 = τ2 3 4τ4 , ϑ = α+γℵ+k(β+kµ)+CΩ2+τ2 Ω . Applying the set (5.1) of solutions, some analytical solutions for Eq. (1) are obtained as follows: (5.1,1) If τ2 < 0 and ϱ ̸= 0, the below singular periodic solution is raised: ϕ5.1,1 = 6τ2 ϱ csc2 [ (x+ ky + ℵz − Ωt) √ −τ2 ] . (40) (5.1,2) If τ2 > 0, τ4 > 0, τ3 ̸= 2ε √ τ2τ4 where ε = ±1 and ϱ ̸= 0, a singular soliton type of solution is obtained as follows: ϕ5.1,2 = −6τ2 ϱ csch2 [(x+ ky + ℵz − Ωt) √ τ2] . (41) By inserting the set of solutions (5.2) in Eq. (1, the following solution can be obtained: (5.2) If τ2 > 0, τ4 > 0, τ3 = −2 √ τ2τ4 and ϱ ̸= 0, the following bright soliton solution: ϕ5.2 = 3τ2 2ϱ sech2 [ (x+ ky + ℵz − Ωt) √ τ2 4 ] . (42) W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 13 of 23 4. Results and Discussion By changing the parameters in the model under investigation, several sets of values that had been documented or reached for Eq. (1) were found. This section includes a variety of graph forms, such as contour plots, 3-D, and 2-D plots of several individual solutions, to help illuminate the mathematical and physical characteristics of the retrieved solutions. Figure (2) shows the dark soliton solution representation of Eq. (12), when selecting τ2 = −0.8, k = 0.7, ℵ = 0.6, Ω = −0.45, α = 0.9, β = 1.9, µ = 0.9, ϑ = 0.6, C = 0.5, γ = 0.5, ϱ = −0.5, y = z = 0, 0 ≤ t ≤ 5 and −30 ≤ x ≤ 30. Localized areas of reduced amplitudes within a surrounding medium are known as dark solitons [30]. In fluid dynamics, they could be connected to regions of decreased fluid density or pressure [31]. The dark soliton solutions, which usually appear in defocusing nonlinear media, are localized intensity dips contained in a continuous wave backdrop. They are important in the context of Bose-Einstein condensates and optical communications, and they simulate phase- shifted energy voids. By settings the values of τ2 = −0.8, k = 0.7, ℵ = 0.6, Ω = −0.45, ϱ = 3, y = z = 0, 0 ≤ t ≤ 10 and −15 ≤ x ≤ 15. This yields the singular periodic solution of Eq. (40) as shown in Figure (3). In the context of physical phenomena, a system is referred to as having a unique periodic solution when there exist both periodic activity and sudden changes or extreme occurrences. These discontinuities or singularities might be brought on by external stimuli, boundary conditions, or non-linearities. Figure (4) displays the bright soliton solution of Eq. (42) with setting the parameters to be as τ2 = 0.8, k = 0.9, ℵ = 0.8, Ω = −0.6, ϱ = 0.9, y = z = 0, 0 ≤ t ≤ 5 and −15 ≤ x ≤ 15. The bright soliton solution has a localized intensity peak above a continuous wave background, thus, it can exist on a finite-depth fluid in the water depth range where the carrier waves have a stable modulation.Because of a careful balancing act between dispersive spreading and nonlinear self-focusing, the resulting dazzling soliton solutions correspond to localized wave packets that maintain their form. These structures are essential for energy localization and distortion-free transmission and are typical of focusing nonlinear media, including optical fibers and certain plasma environments. Figure (5) illustrates the singular soliton solution of Eq. (41) with parameters τ2 = 0.8, k = 0.9, ℵ = 0.8, Ω = −0.6, ϱ = 2.9, y = z = 0, 0 ≤ t ≤ 5 and −15 ≤ x ≤ 15. The NPDE solutions that exhibit singular behavior are known as singular solitons, which are mainly characterized by a narrow region due to an infinity approach peak. They may serve as symbolic representations of local events. Because they are so severe, they are less prevalent in physical systems. A divergence in wave amplitude is shown by the singular soliton solutions’ unique behavior at particular times or locations in space. These solutions frequently simulate wave breaking in very nonlinear dispersive media, energy collapse occurrences, and shock-like structures. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 14 of 23 (a) 3-D (b) contour (c) 2-D Figure 2: Visualization of the solution in Eq. (12). 5. Comparison with literature In Table 1, we present a novel comparative analysis, supported by recent references, from multiple perspectives, including the employed methods, the nature of the obtained solutions, and the specific advantages of each approach when applied to the same class of models investigated in this work. 6. Linear Stability Analysis In this section, we investigate the linear stability of the steady-state solution of the governing equation. 6.1. Steady-State and Perturbed Solution The steady-state solution is one that exhibits no variation in time [36]. For simplicity, we assume it to be a constant, denoted by λ. We then consider a perturbed solution of the form: ϕ(x, y, z, t) = λ+ ρω(x, y, z, t), where ω(x, y, z, t) is a small perturbation function, and ρ is a small parameter. Substituting this form into the original equation (Eq. 1) and neglecting nonlinear terms yields the linearized equation: Cωtt + ωxxxx + αωxx + β ωxy + γ ωxz + ϑωxt + µωyy + 2λϱωxx = 0. (43) W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 15 of 23 (a) 3-D (b) contour (c) 2-D Figure 3: Visualization of the solution in Eq. (40). 6.2. Normal Mode Analysis To analyze Eq. (43), we assume a solution of the form [37]: ω(x, y, z, t) = η ei(Mx+Ry+Fz+Wt), where: • M , R, and F are the wave numbers in the x, y, and z directions, respectively, • W is the temporal frequency, • η is the amplitude. Substituting this ansatz into Eq. (43) gives the algebraic dispersion relation: M4 − CW 2 − αM2 − µR2 −M(γF + βR+ ϑW )− 2M2λϱ = 0. (44) 6.3. Dispersion Relation Solving for W , we obtain the frequency as a function of the wave numbers: W = −ϑM ± √ −4CM(γF + βR) + 4CM4 +M2 (ϑ2 − 4C(α+ 2λϱ))− 4CµR2 2C (45) This relation characterizes the temporal evolution of the perturbations and can be used to assess the stability of the steady-state solution. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 16 of 23 (a) 3-D (b) contour (c) 2-D Figure 4: Visualization of the solution in Eq. (42). 6.4. Compact Stability Criterion Recall the dispersion relation for mode (M,R,F ): W = −ϑM ± √ ∆(M,R,F ) 2C , ∆(M,R,F ) = −4CM(γF+βR)+4CM4+M2 ( ϑ2 − 4C(α+ 2λϱ) ) −4CµR2. Write W = Wre + iWim. The time dependence is eiWt = eiWrete−Wimt. Thus: (i) If Wim > 0, mode decays ⇒ stable. (ii) If Wim < 0, mode grows ⇒ unstable. (iii) If Wim = 0, neutrally stable (pure oscillation). Since for real (M,R,F ), ∆ is real: • Case ∆ ≥ 0: √ ∆ ∈ R ⇒ W ∈ R ⇒ Wim = 0. Modes are neutrally stable. • Case ∆ < 0: In this case, we write ∆ = −|∆|, so the square root becomes imaginary:√ ∆ = i √ |∆|. Then the frequency becomes W = −ϑM 2C ± i √ |∆| 2C . Hence, the imaginary part is Wim = ± √ |∆| 2C . W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 17 of 23 (a) 3-D (b) contour (c) 2-D Figure 5: Visualization of the solution in Eq. (41). Therefore: – If C > 0, the solution to the linearized equation (43) can be expressed as a superposition of two complex conjugate modes corresponding to the roots of the dispersion equation. That is, ω(x, y, z, t) = e i ( Mx+Ry+Fz+ ( −ϑM 2C +i √ |∆| 2C ) t ) + e i ( Mx+Ry+Fz+ ( −ϑM 2C −i √ |∆| 2C ) t ) . Even though one of the branches corresponds to a decaying mode, the other represents an exponentially growing mode. As time evolves, this growing mode will dominate the dynamics, rendering the steady-state solution linearly unstable. Therefore, the instability is invitable. – The existence of even a single wave mode (M,R,F ) for which ∆ < 0 is sufficient to conclude linear instability of the steady state. In the following, we present a set of plots illustrating both instability and neutral stability cases based on different parameter settings. For the neutrally stable case, we adopt the parameter set: α = 1, β = 5, γ = 5, C = 1, F = 0, µ = 1, ϑ = 1, R = 2, ϱ = 5. These values ensure that the discriminant remains non-negative, resulting in a purely real dispersion relation and therefore neutral stability. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 18 of 23 Article Ref. No. Method Outcomes Advantages [32] Bilinear method Multi-soliton solutions, bright and dark solitons, soliton inter- actions Highly effective for integrable systems; systematic and elegant for constructing multi-soliton solutions [33] Sine–cosine method Periodic solutions, bright/dark solitons in trigonometric or hy- perbolic form Simple implementation; effec- tive for constructing wave-type solutions with few parameters [34] Lie symmetry method Similarity reductions, invariant solutions, symmetry classifica- tion Powerful for finding reduction to ODEs, conserved quantities, and hidden symmetries [35] Darboux transfor- mation Bright/dark solitons, rogue and breather solutions Constructs exact solutions re- cursively; excellent for inte- grable models with Lax pairs In this work Improved Modified Extended Tanh- Function Method Bright, dark, and singular soli- tons; periodic, rational, JEFSs, and Weierstrass doubly elliptic solutions Unifies diverse solution types; handles nonlinearities system- atically; suitable for higher- dimensional models Table 1: Comparison of classical methods with the IME tanh-function technique used in this study for (3+1)-dimensional models For the instability scenario, we consider the parameters: α = 3, β = 3, γ = −3, C = 1, F = 1, µ = 2, ϑ = 1, R = 1, ϱ = 1. In this case, the discriminant ∆ < 0, leading to a nonzero imaginary part of the dispersion relation. We therefore plot the imaginary part of ω(M) versus the wavenumber M to identify the regions of instability. 6.5. Conclusion on Stability The steady-state solution λ is: • Linearly unstable if there exists any real mode (M,R,F ) such that ∆(M,R,F ) < 0. In this case, the corresponding perturbation grows exponentially with time for one branch of the root, and the other branch of the root implying decaying, resulting in overall instability of the steady state. • Linearly (neutrally) stable if ∆(M,R,F ) ≥ 0 for all real wave numbers (M,R,F ). Then, all perturbations are purely oscillatory and remain bounded for all time, but they do not decay. This analysis shows that for linear stability (i.e., decay of perturbations), one would require Wim < 0 for all modes, which does not occur in the current conservative system. Therefore, at best, the steady state is marginally (neutrally) stable when ∆ ≥ 0 for all modes. 7. Conclusions This investigation was conducted of the (3+1)-dimensional Boussinesq-KP-type equation, which is commonly employed in models of physical phenomena such as fluid dynamics, physics, hydrody- namic models, and optical mathematical models. This equation is more accurate than either the W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 19 of 23 (a) neutrally stable case (b) instability case Figure 6: Representing the neutrally stable and instability cases for different values of λ KPE or Boussinesq-type equation, each separately, for the traveling waves. A well-known analyti- cal method, the IME tanh-function algorithm, was applied, for which multi-soliton solutions and other solutions were obtained. The resulting solutions included dark, bright, and singular solitons, singular periodic solutions, JEFSs, rational solutions, exponential solutions, moreover Weierstrass elliptic doubly periodic solutions. This study offered a comprehensive explanation for numer- ous intricate and nonlinear phenomena that emerge and spread throughout a variety of nonlinear physical and engineering systems, particularly the nonlinear phenomena found in fluid mechanics, seas, and oceans, not to mention the abundance of nonlinear phenomena found in plasma physics. In contrast to other approaches, the strategy of this study provided new insights into the issues with this model’s problems. This study was potentially expanded to investigate soliton dynamics in multi-dimensional systems, as it became easier to examine through the offered 2D, 3D, and contour graphs for some of the retrieved solutions. Furthermore, a linear stability assessment is conducted to examine the robustness and persistence of the derived solutions under small perturbations. Acknowledgements This research has been funded by Scientific Research Deanship at the University of Ha’il-Saudi Arabia through project number RG-25005. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 20 of 23 References [1] Bang-Qing Li, Abdul-Majid Wazwaz, and Yu-Lan Ma. Two new types of nonlocal boussinesq equations in water waves: bright and dark soliton solutions. Chinese Journal of Physics, 77:1782–1788, 2022. [2] Gui-Qiong Xu and Abdul-Majid Wazwaz. A new (n+ 1)-dimensional generalized kadomtsev– petviashvili equation: integrability characteristics and localized solutions. Nonlinear Dynam- ics, 111(10):9495–9507, 2023. [3] Karim K Ahmed, Niveen M Badra, Hamdy M Ahmed, and Wafaa B Rabie. Unveiling op- tical solitons and other solutions for fourth-order (2+ 1)-dimensional nonlinear schrödinger equation by modified extended direct algebraic method. Journal of Optics, pages 1–13, 2024. [4] SA Khuri and Abdul-Majid Wazwaz. Optical solitons and traveling wave solutions to kudryashov’s equation. Optik, 279:170741, 2023. [5] Muhammad Amin S Murad, Mohammed A Mustafa, Usman Younas, Homan Emadifar, Abeer S Khalifa, Wael W Mohammed, and Karim K Ahmed. Soliton solutions to the gener- alized derivative nonlinear schrödinger equation under the effect of multiplicative white noise and conformable derivative. Scientific Reports, 15(1):1–15, 2025. [6] M Elsaid Ramadan, Hamdy M Ahmed, Abeer S Khalifa, and Karim K Ahmed. Invariant solitons and travelling-wave solutions to a higher-order nonlinear schrödinger equation in an optical fiber with an improved tanh-function algorithm. Journal of Applied Analysis & Computation, 15(6):3270–3289, 2025. [7] Wafaa B Rabie, Karim K Ahmed, Niveen M Badra, Hamdy M Ahmed, M Mirzazadeh, and M Eslami. New solitons and other exact wave solutions for coupled system of perturbed highly dispersive cgle in birefringent fibers with polynomial nonlinearity law. Optical and Quantum Electronics, 56(5):875, 2024. [8] Islam Samir, Hamdy M Ahmed, Homan Emadifar, and Karim K Ahmed. Traveling and soliton waves and their characteristics in the extended (3+ 1)-dimensional kadomtsev–petviashvili equation in fluid. Partial Differential Equations in Applied Mathematics, 14:101146, 2025. [9] Abeer S Khalifa, Hamdy M Ahmed, Niveen M Badra, and Wafaa B Rabie. Exploring solitons in optical twin-core couplers with kerr law of nonlinear refractive index using the modified extended direct algebraic method. Optical and Quantum Electronics, 56(6):1060, 2024. [10] Yinshen Xu, Dumitru Mihalache, and Jingsong He. Resonant collisions among two- dimensional localized waves in the mel’nikov equation. Nonlinear Dynamics, 106:2431–2448, 2021. [11] John Weiss. The painlevé property for partial differential equations. ii: Bäcklund transforma- tion, lax pairs, and the schwarzian derivative. Journal of Mathematical Physics, 24(6):1405– 1413, 1983. [12] Abdul-Majid Wazwaz. Partial differential equations. CRC Press, 2002. [13] Abeer S Khalifa, Hamdy M Ahmed, Niveen M Badra, Jalil Manafian, Khaled H Mahmoud, Kottakkaran Sooppy Nisar, and Wafaa B Rabie. Derivation of some solitary wave solutions for the (3+ 1)-dimensional pkp-bkp equation via the ime tanh function method. AIMS Math- ematics, 9(10):27704–27720, 2024. [14] Ryogo Hirota. The direct method in soliton theory. Number 155. Cambridge university press, 2004. [15] Mohammed Bakheet Almatrafi. Solitary wave solutions to a fractional model using the im- proved modified extended tanh-function method. Fractal and Fractional, 7(3):252, 2023. [16] MB Almatrafi. Construction of closed form soliton solutions to the space-time fractional sym- metric regularized long wave equation using two reliable methods. Fractals, 31(10):2340160, 2023. W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 21 of 23 [17] Muhammad Amin S Murad, Hajar F Ismael, Tukur A Sulaiman, Nehad A Shah, and Jae Dong Chung. Higher-order time-fractional sasa–satsuma equation: Various optical soliton solutions in optical fiber. Results in Physics, 55:107162, 2023. [18] Usman Younas, Tukur A Sulaiman, Hajar F Ismael, and Muhammad Amin S Murad. On the study of interaction phenomena to the (2+ 1)-dimensional korteweg–de vries–sawada–kotera– ramani equation. Modern Physics Letters B, 39(09):2450437, 2025. [19] Usman Younas, Tukur Abdulkadir Sulaiman, and Jingli Ren. Propagation of m-truncated optical pulses in nonlinear optics. Optical and Quantum Electronics, 55(2):102, 2023. [20] Usman Younas, Tukur Abdulkadir Sulaiman, and Jingli Ren. On the optical soliton structures in the magneto electro-elastic circular rod modeled by nonlinear dynamical longitudinal wave equation. Optical and Quantum Electronics, 54(11):688, 2022. [21] Joseph Boussinesq. Essai sur la théorie des eaux courantes. Impr. nationale, 1877. [22] Abdul-Majid Wazwaz. Multiple-soliton solutions for the boussinesq equation. Applied Math- ematics and Computation, 192(2):479–486, 2007. [23] AM Wazwaz. Construction of soliton solutions and periodic solutions of the boussinesq equa- tion by the modified decomposition method. Chaos, Solitons & Fractals, 12(8):1549–1556, 2001. [24] Abdul-Majid Wazwaz. Partial differential equations and solitary waves theory. Springer Science & Business Media, 2010. [25] WEAAM Alhejaili, ABDUL-MAJID Wazwaz, and SA El-Tantawy. On the multiple soliton and lump solutions to the (3+ 1)-dimensional painlev e integrable boussinesq-type and kp-type equations. Rom. Rep. Phys, pages 1–24, 2024. [26] Karim K Ahmed, Niveen M Badra, Hamdy M Ahmed, and Wafaa B Rabie. Soliton solutions and other solutions for kundu–eckhaus equation with quintic nonlinearity and raman effect using the improved modified extended tanh-function method. Mathematics, 10(22):4203, 2022. [27] Mina M Fahim, Hamdy M Ahmed, KA Dib, and Islam Samir. Derivation of dispersive solitons with quadrupled power law of nonlinearity using improved modified extended tanh function method. Journal of Optics, pages 1–10, 2024. [28] Karim K Ahmed, Niveen M Badra, Hamdy M Ahmed, and Wafaa B Rabie. Soliton solutions of generalized kundu-eckhaus equation with an extra-dispersion via improved modified extended tanh-function technique. Optical and Quantum Electronics, 55(4):299, 2023. [29] Zonghang Yang and Benny YC Hon. An improved modified extended tanh-function method. Zeitschrift für Naturforschung A, 61(3-4):103–115, 2006. [30] Hadi Susanto and Magnus Johansson. Discrete dark solitons with multiple holes. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 72(1):016605, 2005. [31] Anne Maitre, Giovanni Lerario, Adrià Medeiros, Ferdinand Claude, Quentin Glorieux, Elisa- beth Giacobino, Simon Pigeon, and Alberto Bramati. Dark-soliton molecules in an exciton- polariton superfluid. Physical Review X, 10(4):041028, 2020. [32] Halide Gümüş and Abdullah Baykal. Exact solutions of boussinesq equations by hirota direct method. Afyon Kocatepe Üniversitesi Fen Ve Mühendislik Bilimleri Dergisi, 24(5):1113–1119, 2024. [33] Sadaf Bibi and Syed Tauseef Mohyud-Din. Traveling wave solutions of kdvs using sine–cosine method. Journal of the Association of Arab Universities for Basic and Applied Sciences, 15(1):90–93, 2014. [34] Sachin Kumar and Shubham Kumar Dhiman. Lie symmetry analysis, optimal system, ex- act solutions and dynamics of solitons of a (3+ 1)-dimensional generalised bkp–boussinesq equation. Pramana, 96(1):31, 2022. [35] Hongyu Luo, Chunxiao Guo, Yanfeng Guo, and Jingyi Cui. Breathing wave solutions and W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 22 of 23 y-type soliton soluions of the new (3+ 1)-dimensional pkp-bkp equation. Nonlinear Dynamics, 112(22):20129–20139, 2024. [36] Mustafa Inc, Abdullahi Yusuf, Aliyu Isa Aliyu, and Dumitru Baleanu. Soliton solutions and stability analysis for some conformable nonlinear partial differential equations in mathematical physics. Optical and Quantum Electronics, 50:1–14, 2018. [37] Mina M Fahim, Hamdy M Ahmed, KA Dib, M Elsaid Ramadan, and Islam Samir. Construct- ing the soliton wave structure and stability analysis to generalized calogero–bogoyavlenskii– schiff equation using improved simple equation method. AIMS Mathematics, 10(5):11052– 11070, 2025. Appendix A In this part, we show all types of solutions for Eq. (6) as shown in [28, 29]. Family (1): τ0 = τ1 = τ3 = 0, W(η) = √ −τ2 τ4 sech( √ τ2 η), τ2 > 0, τ4 < 0, W(η) = √ −τ2 τ4 sec( √ −τ2 η), τ2 < 0, τ4 > 0, W(η) = −ε √ τ4 η , τ2 = 0, τ4 > 0. Family (2): τ1 = τ3 = 0, W(η) = ε √ − τ2 2τ4 tanh (√ −τ2 2 η ) , τ2 < 0, τ4 > 0, τ0 = τ22 4τ4 , W(η) = ε √ τ2 2τ4 tan (√ τ2 2 η ) , τ2 > 0, τ4 > 0, τ0 = τ22 4τ4 , W(η) = √ −τ2m 2 τ4(2m2 − 1) cn (√ τ2 2m2 − 1 η ) , τ2 > 0, τ4 < 0, τ0 = τ22m 2(1−m2) τ4(2m2 − 1)2 , W(η) = √ −m2 τ4(2−m2) dn (√ τ2 2−m2 η ) , τ2 > 0, τ4 < 0, τ0 = τ22 (1−m2) τ4(2−m2)2 , W(η) = ε √ − τ2m 2 τ4(1 +m2) sn (√ − τ2 1 +m2 η ) , τ2 < 0, τ4 > 0, τ0 = τ22m 2 τ4(m2 + 1)2 , where m is the modulus of the Jacobi elliptic functions, 0 ≤ m ≤ 1. Family (3): τ2 = τ4 = 0, τ0 ̸= 0, τ1 ̸= 0, τ3 > 0, A Weierstrass elliptic doubly periodic type solution is obtained: W(η) = ℘ [√ τ3 2 η, A2, A3 ] , W. W. Mohammed et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6490 23 of 23 where A2 = −4f1 f3 and A3 = −4f0 f3 are called Weierstrass elliptic function in-variants. Family (4): τ3 = τ4 = 0, W(η) = − τ1 2τ2 + exp (ε √ τ2 η) , τ2 > 0, τ0 = τ21 4τ2 , W(η) = − τ1 2τ2 + ετ1 2τ2 sin (√ −τ2 η ) , τ0 = 0, τ2 < 0, W(η) = − τ1 2τ2 + ετ1 2τ2 sinh (2 √ τ2 η) , τ0 = 0, τ2 > 0, W(η) = ε √ −τ0 τ2 sin (√ −τ2 η ) , τ1 = 0, τ0 > 0, τ2 < 0, W(η) = ε √ τ0 τ2 sinh ( √ τ2 η) , τ1 = 0, τ0 > 0, τ2 > 0. Family (5): τ0 = τ1 = 0, τ4 > 0, W(η) = − τ2 sec 2 (√ −τ2 2 η ) 2ε √ −τ2τ4 tan (√ −τ2 2 η ) + τ3 , τ2 < 0, W(η) = τ2sech 2 (√ τ2 2 η ) 2ε √ τ2τ4 tanh (√ τ2 2 η ) − τ3 , τ2 > 0, τ3 ̸= 2ε √ τ2τ4, W(η) = 1 2 ε √ τ2 τ4 [ 1 + tanh (√ τ2 2 η )] , τ2 > 0, τ3 = 2ε √ τ2τ4.