On the Extended Simple Equations Method (SEsM) for Obtaining Numerous Exact Solutions to Fractional Partial Differential Equations: A Generalized Algorithm and Several Applications

1. Introduction

Over the past two decades, the use of fractional nonlinear partial differential equations (FNPDEs) has gained considerable popularity in the scientific community as a powerful mathematical framework for modeling and analyzing a wide range of complex phenomena observed in the real world. This modeling approach has been successfully applied in numerous scientific fields, including biology and ecology [1,2,3] , fluid mechanics [4,5,6,7] , finance and economics [8,9,10] and engineering [11,12,13] , etc. A central direction in the study of FNPDEs is finding and analyzing their exact analytical solutions. Analytical solutions provide deeper insight into the qualitative behavior of the modeled systems and serve as criteria for validating numerical simulations. However, obtaining exact solutions of FNPDEs poses significant mathematical challenges due to the combined complexities introduced by nonlinearity and non-integer order of derivatives.

Two main strategies for obtaining exact solutions of FNPDEs are widely used in the literature. The first approach involves applying appropriate fractional transformations to reduce the original FNPDEs into equivalent systems of nonlinear ordinary differential equations (ODEs) of integer order. Once this reduction is achieved, various well established analytical techniques can be used, such as ansatz–based methods [14]–19], group analysis [20,21], perturbation techniques [22,23], and decomposition methods [24,25]. Notable among these are the Hirota method [26], the homogeneous balance method [27], the auxiliary equation method [28], the Jacobi elliptic function expansion method [29], the (G’/G)-expansion method [30], the F–expansion method [15], the tanh–function method [32], the first integral method [33], the exponential function method [34], the simplest equation method [35], the modified method of simplest equation [36,37],etc. We note that the SEsM used in this study has an almost universal character, since most existing methods for obtaining analytical solutions of nonlinear PDEs (whether via integer–order ODEs or specific special functions) are limited to its particular cases, as shown in [38,39].

The second approach relies on standard traveling wave transformations that reduce the original FNPDEs to simpler fractional differential equations, usually fractional ODEs. Analytical solutions to these fractional ODEs are often represented by known special functions such as the Mittag-Leffler function. These functions are also involved in constructing solutions to the original FNPDEs. Within this class of methods, the fractional sub-equation method [40]–43] is particularly prominent. It uses known solutions to fractional versions of Riccati ODEs to construct analytical solutions to the FNPDEs under study.

The extended SEsM presented in this work provides a flexible and unified framework for obtaining analytical solutions to both single FNPDEs and systems of FNPDEs. In its classical form, SEsM [44,45,46] was designed to find exact solutions to NPDEs. According to the original SEsM algorithm the analytical solutions of the studied NPDEs are constructed as complex composite functions, which involve one or more simple composite functions that are power series of solutions to one or more simple equations with the same independent variables. In the last few years, an upgrade of SEsM has been proposed in [47,49,50]. It consists of presenting the analytical solutions of NPDEs (FNPDEs) as composite functions involving the solutions of at least two different simple equations with different independent variables. This extended version of SEsM was initially proposed in [47], where it was applied to the extended integer–order Korteweg–de Vries (KdV) equation, with the obtained exact solutions being combinations of two simple composite functions that involve the solutions of two different simple equations with different independent variables. To find the analytical solutions in [47], all popular types of first-order ODEs, such as variants of ODEs of Riccati and Bernoulli and an ODE of Abel of first order, have been used as simple equations. So far, a similar approach has been proposed only in [48]. It is based on the Kudryashov method, and is applied to the integer–order Boussinesq-like system and the integer–order shallow water wave equation. In this study, however, the proposed methodology is limited to using only two specific ODEs (sub–variants of an ODE of Riccati) as simple equations.

Regarding finding exact solutions of systems of FNDEs, two versions of the extended SEsM adapted to systems of FNDEs have recently been proposed in [49,50]. In [49] an extended version of SEsM was proposed based on the assumption that variables in a spatio–temporal dynamical system modeled by FNPDEs propagate at different velocities. In this framework, exact solutions of the studied system are constructed using simple composite functions formed by solutions of similar or different types of simple equations with distinct independent variables. This version of the extended SEsM was applied to a system of two FNPDEs that models a complex ecological phenomenon. In this way, numerous new exact solutions were derived using different solutions of a set of first-order ODEs. In [50] a different version of extended SEsM was introduced under the assumption that variables in the system of FNPDEs exhibit synchronized multi–wave behaviour. Accordingly, the solutions of the studied system are constructed as complex composite functions involving power series of solutions of two or more different types of simple equations with different independent variables. This version of the extended SEsM was applied to a system of two FNPDEs that models a complex fluid–like process, such as the one considered in the current study. In this way, numerous multi–wave solutions were obtained based on known exact solutions of a set of different second-order ODEs.

Building on both versions given above, in this paper, a generalized algorithm of the extended SEsM for obtaining exact solutions to systems of FNPDEs is proposed. The algorithm offers several variants at each step, allowing the user to select the most appropriate variant depending on the expected wave dynamics in the modeled real-world process. In this study, we demonstrate the effectiveness of this algorithm by its application to time-fractional versions of the models presented in [51]:

$$\begin{array}{ccc}& & D_t^{\varphi }{u}_1+{v_1}_x+u_1{u_1}_x-3\alpha {u_1}_{xxt}=0\hfill \\ & & D_t^{\varphi }{v}_1+{u_1}_x+{\left(u_1{v}_1\right)}_x-\alpha {u_1}_{xxx}=0\hfill \end{array}$$

(1)

and

$$\begin{array}{ccc}& & D_t^{\varphi }{u}_2+{v_2}_x+u_2{u_2}_x=0\hfill \\ & & D_t^{\varphi }{v}_2+{\left(u_2{v}_2\right)}_x+{u_2}_{xxx}=0\hfill \end{array}$$

(2)

Equation (1) and Equation (2) are mathematical generalizations inspired by two variants of classical Boussinesq system, where $\varphi$ is a time–fractional number and $\alpha$ in Equation (1) is an arbitrary coefficient. Bellow we shall derive exact solutions of Equation (1) and Equation (2) applying two different step’s variants of the generalized SEsM algoritm regarding the type of the wave transformation introduced and the type of simple equations used. We note that various new exact solutions of a model system similar to Equation (2) were obtained in [50]. However, to derive these analytical solutions, in [50], a different sub–variant of the extended SEsM was applied than the one used in the current study.

2. Preliminaries Relevant to the Present Study

For the convenience of the reader, in this section we present several important statements from [52,53], which clarify some basic aspects in the current study.

2.1. Fractional Derivative via Fractional Difference

Definition 2.1. Let $f:R\to R,x\to f\left(x\right)$, denote a continuous (but not necessarily differentiable) function, and let $h>0$ denote a constant discretization span. Define the forward operator $FW\left(h\right)$ by the equality (the symbol $:=$ means that the left side is defined by the right side)

$$FW\left(h\right)f\left(x\right):=f(x+h);$$

(3)

then the fractional difference of order $\varphi$, $0<\varphi <1$, of $f\left(x\right)$ is defined by the expression

$${\Delta }^{\varphi }f\left(x\right):={(FW-1)}^{\varphi }=\sum _{k=0}^{\infty }{(-1)}^k\left({\displaystyle \genfrac{}{}{0pt}{}{\varphi }{k}}\right)f[x+(\varphi -k)h]$$

(4)

and its fractional derivative of order $\varphi$ is defined by the limit

$$f^{\varphi }\left(x\right)=\underset{h\downarrow 0}{lim}{\displaystyle \frac{{\Delta }^{\varphi }f\left(x\right)}{h^{\varphi }}}$$

(5)

2.2. Modified Fractional Riemann–Liouville Derivative

Definition 2.2. Refer to the function $f\left(x\right)$ of Definition 2.1.

(i) Assume that $f\left(x\right)$ is a constant K . Then its fractional derivative of order $\varphi$ is

$$\begin{array}{c}\hfill D_x^{\varphi }=K{\Gamma }^{-1}(1-\varphi )x^{-\varphi },\varphi \le 0,\\ \hfill =0,\varphi \ge 0.\end{array}$$

(6)

(ii) When $f\left(x\right)$ is not a constant, then we will set

$$f\left(x\right)=f\left(0\right)+\left[f\right(x)-f(0\left)\right],$$

(7)

and its fractional derivative will be defined by the expression

$$f^{\varphi }\left(x\right)=D_x^{\varphi }\left(0\right)+D_x^{\varphi }\left(f\left(x\right)-f\left(0\right)\right),$$

(8)

in which, for negative $\varphi$, one has

$$D_x^{\varphi }(f\left(x\right)-f\left(0\right)):={\displaystyle \frac{1}{{\Gamma }^{-\varphi }}}{\int }_0^x{(x-\zeta )}^{-\varphi -1}f\left(\zeta \right)d\zeta ,\varphi \le 0$$

(9)

whilst for positive $\varphi$ we will set

$$D_x^{\varphi }(f\left(x\right)-f\left(0\right))=D_x^{\varphi }\left(f\left(x\right)\right)=D_x^{\varphi }\left(f^{(\varphi -1)}\left(x\right)\right),\varphi \ge 0.$$

(10)

When $n\le \varphi <n+1$,

$$f^{\varphi }\left(x\right):=f^{(\varphi -1)}{\left(x\right)}^n,n\le \varphi ,n+1$$

(11)

In order to find the fractional derivative of compound functions, the following equation holds:

$$d^{\varphi }f\approx \Gamma (1+\varphi )df,0<\varphi <1.$$

(12)

2.3. Taylor’s Series of Fractional Order

Proposition: Assume that he continuous function $f:R\to R,x\to f\left(x\right)$ has a fractional derivative of order $k\varphi$. For any positive integer k and for any $\varphi$, $0<\varphi <1$, then the following equality holds

$$f(x+h)=\sum _{k=0}^{\infty }{\displaystyle \frac{h^{\varphi k}}{\left(\varphi k\right)!}}{f}^{\left(\varphi k\right)}\left(x\right),$$

(13)

with the notation

$$\Gamma (1+k\varphi )=\left(\varphi k\right)!,0<\varphi \le 1$$

(14)

where $\Gamma (.)$ is the Euler gamma function. Formally, this series can be written

$$f(x+h)=E_{\varphi }\left(h^{\varphi }{D}_x^{\varphi }\right)f\left(x\right)$$

(15)

where $D_x$ is the derivative operator with respect to x, and $E_{\varphi }\left(u\right)$ is the Mittag-Leffler function defined by the expression

$$E_{\varphi }\left(u\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{u^k}{\left(\varphi k\right)!}}=\sum _{k=0}^{\infty }{\displaystyle \frac{u^k}{\Gamma (1+k\varphi )}}$$

(16)

Remark: This fractional Taylor’s series does not hold with the standard Riemann–Liouville derivative, but it is applied to non-differentiable functions only.

Corollary: Assume that $m<\varphi \le m+1,m\to N-0$ and that $f\left(x\right)$ has derivatives of order k (integer), $1\le k\le m$. Assume further that $f^m\left(x\right)$ has a fractional Taylor’s series of order $\varphi -m=:\psi$

provided by the expression

$$f^m(x+h)=\sum _{k=0}^{\infty }{\displaystyle \frac{h^{k(\varphi -m)}}{\Gamma [1+k(\varphi -m\left)\right]}}{D}^{k(\varphi -m)}{f}^{\left(m\right)}\left(x\right),m<\varphi <m+1$$

(17)

Then, integrating this series with respect to h provides

$$f(x+h)=\sum _{k=0}^m{\displaystyle \frac{h^k}{k!}}{f}^{\left(k\right)}\left(x\right)+\sum _{k=0}^{\infty }{\displaystyle \frac{h^{(k\psi +m)}}{\Gamma (k\psi +m+1)}}{f}^{(k\psi +m)}\left(x\right),\psi :=\varphi -m$$

(18)

3. Description of the Generalized Algorithm of the Extended Simple Equations Method (SEsM)

In this section, we shall provide a detailed presentation of the generalized SEsM algorithm, highlighting the various capabilities it offers for obtaining different types of complex analytical solutions of systems of FNPDEs. Several variants of this algorithm can be found in previous works [49,50], and we unite and complement them in this article. As it is shown below, unlike the original SEsM algorithm outlined in [44]–46], we have swapped its first and second steps, incorporating methodological extensions into both steps of the original SEsM algorithm.

Taking into account the problem we set out to solve in this study, we shall demonstrate all the possible steps of the algorithm above mentioned for a system of two FNPDEs with two independent variables:

$$\begin{array}{ccc}& & \Psi (u(x,t),D_t^{\varphi }u(x,t),D_x^{\varphi }u(x,t),D_{tt}^{2\varphi }u(x,t),D_{xx}^{2\varphi }u(x,t),...)=0\hfill \\ & & \Omega (v(x,t),D_t^{\varphi }v(x,t),D_x^{\varphi }v(x,t),D_{tt}^{2\varphi }v(x,t),D_{xx}^{2\varphi }v(x,t),...)=0,0<\varphi \le 1\hfill \end{array}$$

(19)

where D denotes the arbitrary fractional order derivative operator, as superscripts $\varphi ,2\varphi ..$ give its fractional number and subscripts denote time and partial derivatives. Here, $\Psi$ and $\Omega$ are polynomials of $u(x,t)$ and $v(x,t)$ and its derivatives, respectively, where $u=u(x,t)$ and $v=v(x,t)$ are unknown functions.

The generalized extended SEsM algorithm which is adopted for finding exact solutions of FNPDEs includes the following basic steps:

• Construction of the solution of Equation (19). In addition to the known conventional methods for obtaining exact solutions of Equation (19), where such solutions are constructed by power series of the solutions of one simple (auxiliary) equation (or one special function) with the same independent variable for the both system variables in (19), the SEsM provides also several alternative variants:
  • Variant 1: Constructing the solution of Equation (19) by single composite functions. These single composite functions can be:

    (a)

    with distinct independent variables: This solution variant is applicable to real–world dynamical models of a type (19), where it is expected that the system variables demonstrate different wave behaviour and they move with different wave speeds. Thus, a such solution takes the form:

    $$\begin{array}{c}\hfill u\left({\xi }_1(x,t)\right)=F_1\left({\xi }_1(x,t)\right),\\ \hfill v\left({\xi }_2(x,t)\right)=F_2\left({\xi }_2(x,t)\right),\end{array}$$

    (20)

    where ${\xi }_1(x,t),{\xi }_2(x,t)$ are generalized wave variables (to be defined in the next step) and $F_1,F_2$ are simple composite functions of the form:

    $$F_1\left({\xi }_1\right)=\sum _{i_1=1}^{n_1}{ϵ}_{i_1}{\left[f_1\left({\xi }_1\right)\right]}^{i_1},F_2\left({\xi }_2\right)=\sum _{i_2=1}^{n_2}{\eta }_{i_2}{\left[f_2\left({\xi }_2\right)\right]}^{i_2},$$

    (21)

    or

    $$F_1\left({\xi }_1\right)=\sum _{i_1=-n_1}^{n_1}{ϵ}_{i_1}{\left[f_1\left({\xi }_1\right)\right]}^{i_1},F_2\left({\xi }_2\right)=\sum _{i_2=-n_2}^{n_2}{\eta }_{i_2}{\left[f_2\left({\xi }_2\right)\right]}^{i_2},$$

    (22)

    where $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ being solutions of simple equations (to be selected at a later stage of the algorithm), and ${ϵ}_{i_1},{\eta }_{i_2}$ are coefficients to be determined later.

    Remark: The simple equations used may have the same analytical form as shown in [49], but they may have a different form. The form of the simple equations used depends on the physical wave characteristics of the real-world system being modeled.

    (b)

    with a single common wave variable: This solution variant is applicable to real–world dynamical models of a type (19), where it is expected that the system variables demonstrate different wave behaviour but it is synchronized. For this case, the solution of Equation (19) reduces to:

    $$\begin{array}{c}\hfill u\left(\xi (x,t)\right)=F_1\left(\xi (x,t)\right),\\ \hfill v\left(\xi (x,t)\right)=F_2\left(\xi (x,t)\right),\end{array}$$

    (23)

    where $F_1,F_2$ are again two simple composite functions of the form similar to (21) or (22), that include the solutions of two distinct simple equations but with common wave variable.

    Remark: When the simple equations used have an identical analytical form, the solution (23) has the same form as that used in all other known methods for finding exact solutions to systems of FNPDEs (NPDES) to date.

  • Variant 2: Constructing the solution of Equation (19) by complex composite functions including combinations of at least two single composite functions. The single composite functions can be:

    (a)

    with distinct independent variables: This solution variant is applicable to real–world dynamical models of a type (19), where it is expected that the both system variables can exhibit both a synchronized multi–wave behavior or non–synchronized multi–wave behavior as the different waves move with different speeds. The simplest examples of synchronized multi–wave behaviour of the variables of system (19) can be expressed analytically as:

    $$\begin{array}{c}\hfill u(x,t)=F_1\left({\xi }_1(x,t)\right)+F_2\left({\xi }_2(x,t)\right),\\ \hfill v(x,t)=F_3\left({\xi }_1(x,t)\right)+F_4\left({\xi }_2(x,t)\right)\end{array}$$

    (24)

    $$\begin{array}{c}\hfill u(x,t)=F_1\left({\xi }_1(x,t)\right)F_2\left({\xi }_1(x,t)\right),\\ \hfill v(x,t)=F_3\left({\xi }_1(x,t)\right)F_4\left({\xi }_2(x,t)\right),\end{array}$$

    (25)

    $$\begin{array}{c}\hfill u(x,t)={\displaystyle \frac{F_1\left({\xi }_1(x,t)\right)}{F_2\left({\xi }_2(x,t)\right)}},\\ \hfill v(x,t)={\displaystyle \frac{F_3\left({\xi }_1(x,t)\right)}{F_4\left({\xi }_2(x,t)\right)}},\end{array}$$

    (26)

    where ${\xi }_1(x,t),{\xi }_2(x,t)$ are again generalized wave variables, and the single composite functions $F_1\left({\xi }_1\right),F_2\left({\xi }_2\right),F_3\left(x{i}_1\right),F_4\left({\xi }_2\right)$ are expressed as follows:

    $$\begin{array}{c}\hfill F_1\left({\xi }_1\right)=\sum _{i_1=1}^{n_1}{ϵ}_{i_1}{\left[f_1\left({\xi }_1\right)\right]}^{i_1},F_2\left({\xi }_2\right)=\sum _{i_2=1}^{n_2}{\eta }_{i_2}{\left[f_2\left({\xi }_2\right)\right]}^{i_2},\\ \hfill F_3\left({\xi }_1\right)=\sum _{i_3=1}^{n_3}{\theta }_{i_3}{\left[f_1\left({\xi }_1\right)\right]}^{i_3},F_4\left({\xi }_2\right)=\sum _{i_4=1}^{n_4}{\iota }_{i_4}{\left[f_2\left({\xi }_2\right)\right]}^{i_4},\end{array}$$

    (27)

    where $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ are solutions of two distinct simple equations, that will be defined in Step 3 of the SEsM algorithm.

    In the case of different multi–wave dynamics implemented by the variables of Equation (19), its solutions can be constructed as combinations of the solutions (24–26). For example, several alternative solutions can been written:

    $$\begin{array}{c}\hfill u(x,t)=F_1\left({\xi }_1(x,t)\right)+F_2\left({\xi }_2(x,t)\right),\\ \hfill v(x,t)=F_3\left({\xi }_1(x,t)\right)F_4\left({\xi }_2(x,t)\right),\end{array}$$

    (28)

    $$\begin{array}{c}\hfill u(x,t)=F_1\left({\xi }_1(x,t)\right)F_2\left({\xi }_2(x,t)\right),\\ \hfill v(x,t)=F_3\left({\xi }_1(x,t)\right)+F_4\left({\xi }_2(x,t)\right),\end{array}$$

    (29)

    $$\begin{array}{c}\hfill u(x,t)=F_1\left({\xi }_1(x,t)\right)+F_2\left({\xi }_2(x,t)\right),\\ \hfill v(x,t)={\displaystyle \frac{F_3\left({\xi }_1(x,t)\right)}{F_4\left({\xi }_2(x,t)\right)}},\end{array}$$

    (30)

    $$\begin{array}{c}\hfill u(x,t)={\displaystyle \frac{F_1\left({\xi }_1(x,t)\right)}{F_2\left({\xi }_2(x,t)\right)}},\\ \hfill v(x,t)=F_3\left({\xi }_1(x,t)\right)+F_4\left({\xi }_2(x,t)\right),\end{array}$$

    (31)

    $$\begin{array}{c}\hfill u(x,t)=F_1\left({\xi }_1(x,t)\right)F_2\left({\xi }_2(x,t)\right),\\ \hfill v(x,t)={\displaystyle \frac{F_3\left({\xi }_1(x,t)\right)}{F_4\left({\xi }_2(x,t)\right)}},\end{array}$$

    (32)

    $$\begin{array}{c}\hfill u(x,t)={\displaystyle \frac{F_1\left({\xi }_1(x,t)\right)}{F_2\left({\xi }_2(x,t)\right)}},\\ \hfill v(x,t)=F_3\left({\xi }_1(x,t)\right)F_4\left({\xi }_2(x,t)\right),\end{array}$$

    (33)

    where the single composite functions $F_1\left({\xi }_1\right),F_2\left({\xi }_2\right),F_3\left({\xi }_1\right),F_4\left({\xi }_2\right)$ are given in Equation (27).

    (b)

    with a single common wave variable: This solution variant is applicable to real–world dynamical models of a type (19), where it is expected that the system variables demonstrate multi–wave behaviour, where the waves propagate with the same speed. In this scenario, ${\xi }_1={\xi }_2=\xi$. Given this, it is easy to make a change in the wave coordinates in Equations (24)–(33) to obtain analogous variants of analytical solutions of Equation (19) for this specific case. The form of some solution variants in this category approaches the forms of analytical solutions proposed in several similar methodologies in this field.

    Remark: The analytical forms of the single composite functions can be of the same type for the both system variables, i.e. they can include combinations of solutions of simple equations with an identical analytical form. However, they can also include combinations of solutions of different types’ simple equations, as the specific construction forms being determined by the specific physical nature of the model equations.

• Selection of the traveling wave type transformation. To apply the SEsM to Equation (19), it is crucial to define the fractional derivatives in those equations. The choice of fractional derivatives (e.g., Riemann-Liouville, Caputo,conformable, etc.) is essential for accurately modeling wave dynamics and reflecting the system’s physical properties, based on factors like the process nature, boundary conditions, and memory effects interpretation. In this context, the following variants of transformations are possible:
  • Variant 1: Use of a fractional transformation. The choice of explicit form of the fractional traveling wave transformation depends on how the fractional derivatives in Equation (19) are defined. Bellow, the most used fractional traveling wave transformations are selected:

    (a)

    Conformable fractional traveling wave transformation: ${\xi }_1={\kappa }_1{\displaystyle \frac{x^{\varphi }}{\varphi }}+{\omega }_1{\displaystyle \frac{t^{\varphi }}{\varphi }},{\xi }_2={\kappa }_2{\displaystyle \frac{x^{\varphi }}{\varphi }}+{\omega }_2{\displaystyle \frac{t^{\varphi }}{\varphi }}$, defined for conformable fractional derivatives [54];

    (b)

    Fractional complex transform: ${\xi }_1={\kappa }_1{\displaystyle \frac{x^{\varphi }}{\Gamma (\varphi +1)}}+{\omega }_1{\displaystyle \frac{t^{\varphi }}{\Gamma (\varphi +1)}},{\xi }_2={\kappa }_2{\displaystyle \frac{x^{\varphi }}{\Gamma (\varphi +1)}}+{\omega }_2{\displaystyle \frac{t^{\varphi }}{\Gamma (\varphi +1)}}$, defined for modified Riemann–Liouville fractional derivatives [52], which can applied for Caputo fractional derivatives and other fractional derivative types in studied FNPDEs [55].

    (c)

    Some more specific fractional transformations, such as ${\xi }_1={\displaystyle \frac{{\kappa }_1}{\varphi }}{(x+{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}})}^{\varphi }+{\displaystyle \frac{{\omega }_1}{\varphi }}{(t+{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}})}^{\varphi },{\xi }_2={\displaystyle \frac{{\kappa }_2}{\varphi }}{(x+{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}})}^{\varphi }+{\displaystyle \frac{{\omega }_2}{\varphi }}(t+{\displaystyle \frac{1}{\Gamma {\left(\varphi \right)}^{\varphi }}}$, defined for Beta derivatives [56]; ${\xi }_1={\displaystyle \frac{\Gamma (\Theta +1)}{\varphi }}({\kappa }_1{x}^{\varphi }+{\omega }_1{t}^{\varphi }),{\xi }_2={\displaystyle \frac{\Gamma (\Theta +1)}{\varphi }}({\kappa }_2{x}^{\varphi }+{\omega }_2{t}^{\varphi })$, defined for M-truncated derivatives [57] and others.

    In all the cases, the studied FNPDEs are reduced to integer–order nonlinear ODEs.
  • Variant 2: Use a standard traveling wave transformation. In this case, by introducing a traveling wave ansatz ${\xi }_1={\kappa }_{1}x+{\omega }_{1}t,{\xi }_2={\kappa }_{2}x+{\omega }_{2}t$ in the selected variant solutions from Step 1, the studied FNPDEs are reduced to fractional nonlinear ODEs.
• Selection of the forms of the used simple equations.
  • For Variant 1 of Step 2: The general form of the integer–order simple equations used is:

    $${\left({\displaystyle \frac{d^{k_1}{f}_1}{d{\xi }_1^{k_1}}}\right)}^{l_1}=\sum _{j_1=0}^{m_1}{\mu }_{j_1}{f}_1^{j_1},{\left({\displaystyle \frac{d^{k_2}{f}_2}{d{\xi }_2^{k_2}}}\right)}^{l_2}=\sum _{j_2=0}^{m_2}{\nu }_{j_2}{f}_2^{j_2}.$$

    (34)

    By fixing $l_1,l_2$ and $m_1,m_2$ (at $k_1=k_2=1$) in Equation (34), different types integer–order ODEs can be used as simple equations, such as:

    -

    (a) ODEs of first order with known analytical solutions (For example, an ODE of Riccati, an ODE of Bernoulli, an ODE of Abel of first kind, an ODE of tanh–function,etc);

    -

    (b) ODEs of second order with known analytical solutions (For example elliptic equations of Jaccobi and Weiershtrass and their sub–variants, an ODE of Abel of second kind, etc.).

  • For Variant 2 of Step 2: The general form of the fractional simple equations used is:

    $$D_{{\xi }_1}^{\varphi }=\sum _{j_1=0}^{m_1}{\mu }_{j_1}{f}_1^{j_1},D_{{\xi }_2}^{\varphi }=\sum _{j_2=0}^{m_2}{\nu }_{j_2}{f}_2^{j_2},$$

    (35)

    By fixing and $m_1,m_2$ in Equation (35), several types fractional ODEs can be used as simple equations, such as a fractional ODE of Riccati (and its sub–variants) and a fractional ODE of Bernoulli.
• Derivation of the balance equations and the system of algebraic equations. The fixation of the explicit form of constructed variant solutions of Equation (19) presented in Step 1 of the SEsM algorithm depends on the balance equations derived. Substitutions of the selected variants from Steps 1,2 and 3 in Equation (19) leads to obtaining polynomials of the functions $f_1$ and $f_2$. The coefficients in front of these functions include the coefficients of the solution of the considered FNPDEs as well as the coefficients of the simple equations used. Analytical solutions of Equation (19) can be extracted only if each coefficient in front of the functions $f_1$ and $f_2$ contains almost two terms. Equating these coefficients to zero leads to formation of a system of nonlinear algebraic equations for each variant chosen according Steps 1,2 and 3 of the SEsM algorithm.
• Derivation of the analytical solutions. Any non–trivial solution of the algebraic system above mentioned leads to a solution of the studied FNPDEs by replacing the specific coefficients in the corresponding variant solutions, given in Step 1 as well as by changing the traveling wave coordinates chosen by the variants given in Step 2. For simplicity, these solutions are expressed through special functions. For a Variant 1 of the Step 3, these special functions are $V_{{\mu }_0,{\mu }_1,...,{\mu }_{j_1}}({\xi }_1;k_1,l_1,m_1)$ and $V_{{\nu }_0,{\nu }_1,...,{\nu }_{j_2}}({\xi }_2;k_2,l_2,m_2)$, as their explicit forms are determined on the basis of the specific form of the simple equations chosen (For reference, see Equation (34)). For a Variant 2 of the Step 3, the special functions are $V_{{\mu }_0,{\mu }_1,...,{\mu }_{j_1}}({\xi }_1,\varphi ;m_1)$ and $V_{{\nu }_0,{\nu }_1,...,{\nu }_{j_2}}({\xi }_2,\varphi ;m_2)$, whose exact forms are determined by the type of fractional simple equations used (For reference, see see Equation (35)).

The generalized algorithm presented above is designed to support researchers in this field to extract analytical solutions of models of type (19), which describe dynamics of natural processes. Using this algorithm, researchers can select among the variants provided in Steps 1, 2, and 3, making their choice on the basis of the specific characteristics of the wave dynamics which is expected for the natural phenomenon being modeled.

In the remainder of this article, we apply two different variants of the generalized SEsM algorithm described above to Equation (1) and Equation (2).

4. Derivation of Multi–Wave Exact Solutions of Equation (1) Using a Fractional Wave Transformation

Following the generalized SEsM algorithm, presented above, we choose to present the solution of Equation (1) by complex composite functions, which include sums of two single simple composite functions with different independent variables, i.e.

$$\begin{array}{ccc}& & u_1({\xi }_1,{\xi }_2)=\sum _{i_1=0}^{n_1}{ϵ}_{i_1}{\left[f_1\left({\xi }_1\right)\right]}^{i_1}+\sum _{i_2=0}^{n_2}{\eta }_{i_2}{\left[f_2\left({\xi }_2\right)\right]}^{i_2}\hfill \\ & & v_1({\xi }_1,{\xi }_2)=\sum _{i_3=0}^{n_3}{\theta }_{i_3}{\left[f_1\left({\xi }_1\right)\right]}^{i_3}+\sum _{i_4=0}^{n_4}{\iota }_{i_4}{\left[f_2\left({\xi }_2\right)\right]}^{i_4}\hfill \end{array}$$

(36)

We choose to define the system (1) by modified Riemann–Liouville fractional derivatives. Then, the fractional wave transformations are:

$${\xi }_1=x+{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (\varphi +1)}},{\xi }_2=x+{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (\varphi +1)}}$$

(37)

We choose to use ODEs of second order as simple equations. We fix their polynomial parts up to fourth degree, i.e.

$${\left({\displaystyle \frac{d{f}_1}{d{\xi }_1}}\right)}^2=\sum _{j_1=0}^4{\mu }_{j_1}{\left[f_1\left({\xi }_1\right)\right]}^{j_1},{\left({\displaystyle \frac{d{f}_2}{d{\xi }_2}}\right)}^2=\sum _{j_2=0}^4{\nu }_{j_2}{\left[f_2\left({\xi }_2\right)\right]}^{j_2}$$

(38)

There are many sub–variants of Equation (38) depending on the numerical values of coefficients ${\mu }_0,..,{\mu }_4$ and ${\nu }_0,..,{\nu }_4$. In this study, we chose to fix ${\mu }_1={\mu }_3=0$ and ${\nu }_0={\nu }_1={\nu }_3=0$. Then, the simple equations used are reduced to

$${\left({\displaystyle \frac{d{f}_1}{d{\xi }_1}}\right)}^2={\mu }_0+{\mu }_2{f}_1^2+{\mu }_4{f}_1^4,{\left({\displaystyle \frac{d{f}_2}{d{\xi }_2}}\right)}^2={\nu }_2{f}_2^2+{\nu }_4{f}_2^4,$$

(39)

and their analytical solutions (used in this study), are presented in Appendix A.1

A balance procedure is performed, and the resulting balance equations are $n_1=n_3=2m_1-2$ and $n_2=n_4=2m_2-2$ Based on this, the general solution of Equation (1), presented by the special functions $V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}({\xi }_1,1,2,4)$ and $V_{0,0,{\nu }_2,0,{\nu }_4}({\xi }_2,1,2,4)$ is

$$\begin{array}{ccc}& & u_1({\xi }_1,{\xi }_2)={ϵ}_0+{ϵ}_1{V}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}({\xi }_1;1,2,4)+{ϵ}_2{V}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0+{\eta }_1{V}_{0,0,{\nu }_2,0,{\nu }_4}({\xi }_2;1,2,4)+{\eta }_2{V}_{0,0,{\nu }_2,0,{\nu }_4}^2({\xi }_2;1,2,4)\hfill \\ & & v_1({\xi }_1,{\xi }_2)={\theta }_0+{\theta }_1{V}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}({\xi }_1;1,2,4)+{\theta }_2{V}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\iota }_0+{\iota }_1{V}_{0,0,{\nu }_2,0,{\nu }_4}({\xi }_2;1,2,4)+{\eta }_2{V}_{0,0,{\nu }_2,0,{\nu }_4}^2({\xi }_2;1,2,4)\hfill \end{array}$$

(40)

Substitution of Equation (36), Equation (37) and Equation (39) in Equation (1) leads to the following system of nonlinear algebraic equations:

$$\begin{array}{cc}\hfill & 4{{ϵ}_2}^2-120\alpha {ϵ}_2{\mu }_4-8{ϵ}_2{\theta }_2=0\hfill \\ \hfill & -30\alpha {ϵ}_1{\mu }_4+6{ϵ}_1{ϵ}_2-6{ϵ}_2{\theta }_1-6{ϵ}_1{\theta }_2=0\hfill \\ \hfill & 4{ϵ}_2{\eta }_2-4{\theta }_2{\eta }_2-4{ϵ}_2{\iota }_2=0\hfill \\ \hfill & -2{ϵ}_2{\iota }_1+2{ϵ}_2{\eta }_1-2{\theta }_2{\eta }_1=0\hfill \\ \hfill & -4{ϵ}_2{\iota }_1-4{\theta }_2{\eta }_1+4{\eta }_1{ϵ}_2=0\hfill \\ \hfill & -2{ϵ}_2+4{\theta }_2+2{{ϵ}_1}^2+2{\omega }_1{ϵ}_2+2{\omega }_2{ϵ}_2+4{ϵ}_0{ϵ}_2+4{\eta }_0{ϵ}_2\hfill \\ \hfill & -2{\omega }_1{\theta }_2-2{\omega }_2{\theta }_2-4{ϵ}_0{\theta }_2-4{ϵ}_1{\theta }_1-4{\eta }_0{\theta }_2-4{\theta }_0{ϵ}_2-4{\iota }_0{ϵ}_2-40\alpha {ϵ}_2{\mu }_2=0\hfill \\ \hfill & 2{ϵ}_1{\eta }_1-2{ϵ}_1{\iota }_1-2{\theta }_1{\eta }_1=0\hfill \\ \hfill & -4{ϵ}_1{\iota }_2-4{\theta }_1{\eta }_2=0\hfill \\ \hfill & 2{\eta }_2{ϵ}_1-2{ϵ}_1{\iota }_2-2{\theta }_1{\eta }_2=0\hfill \\ \hfill & -{ϵ}_1+2{\theta }_1-2{ϵ}_0{\theta }_1-2{\theta }_0{ϵ}_1-5\alpha {ϵ}_1{\mu }_2-2{\iota }_0{ϵ}_1\hfill \\ \hfill & +{\omega }_2{ϵ}_1+2{ϵ}_0{ϵ}_1+{\omega }_1{ϵ}_1-{\omega }_1{\theta }_1-{\omega }_2{\theta }_1+2{\eta }_0{ϵ}_1-2{\eta }_0{\theta }_1=0\hfill \\ \hfill & -120\alpha {\eta }_2{\nu }_4-8{\eta }_2{\iota }_2+4{{\eta }_2}^2=0\hfill \\ \hfill & -6{\eta }_1{\iota }_2+6{\eta }_1{\eta }_2-6{\eta }_2{\iota }_1-30\alpha {\eta }_1{\nu }_4=0\hfill \\ \hfill & 4{ϵ}_0{\eta }_2-2{\omega }_1{\iota }_2+4{\eta }_0{\eta }_2-40\alpha {\eta }_2{\nu }_2+2{\omega }_2{\eta }_2+2{\omega }_1{\eta }_2-4{\eta }_1{\iota }_1-2{\eta }_2\hfill \\ \hfill & +4{\iota }_2-4{ϵ}_0{\iota }_2-4{\theta }_0{\eta }_2-4{\iota }_0{\eta }_2+2{{\eta }_1}^2-2{\omega }_2{\iota }_2-4{\eta }_0{\iota }_2=0\hfill \\ \hfill & -{\eta }_1+2{\iota }_1-2{ϵ}_0{\iota }_1-2{\eta }_0{\iota }_1-5\alpha {\eta }_1{\nu }_2-{\omega }_1{\iota }_1-{\omega }_2{\iota }_1\hfill \\ \hfill & +2{\eta }_0{\eta }_1-2{\theta }_0{\eta }_1-2{\iota }_0{\eta }_1+2{ϵ}_0{\eta }_1+{\omega }_1{\eta }_1+{\omega }_2{\eta }_1=0\hfill \end{array}$$

(41)

One simple non–trivial solution of this system is:

$$\begin{array}{ccc}& & {ϵ}_1={\eta }_1={\theta }_1={\iota }_1=0;{ϵ}_2={\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right);\hfill \\ & & {\eta }_2=-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right);{\theta }_0=0;\hfill \\ & & {\theta }_2=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right);\hfill \\ & & {\iota }_0={\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0;\hfill \\ & & {\iota }_2=-{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right)\hfill \end{array}$$

(42)

where ${ϵ}_0,{\eta }_0,\alpha ,{\mu }_2,{\mu }_4,{\nu }_2,{\nu }_4,{\omega }_1,{\omega }_2$ are free parameters. Substituting the coefficients of Equation (42) in Equation (40) and combining the different analytical solution of Equation (39) leads to the following families of exact solutions of Equation (1):

• Family 1: A solution combined distinct generalized hyperbolic functions with different independent variables:

  $$\begin{array}{ccc}& & u_1^{\left(1\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(1\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (43)

  where

  $$V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(1\right)}({\xi }_1;1,2,4)=\sqrt{{\displaystyle \frac{-{\mu }_2}{2{\mu }_4}}}tanh\left({\displaystyle \frac{{\mu }_2}{2}}{\xi }_1\right)$$

  (44)

  for ${\mu }_0={\displaystyle \frac{{\mu }_2^2}{4{\mu }_4}}$, ${\mu }_2<0,{\mu }_4>0$ and

  $$V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(1\right)}({\xi }_2;1,2,4)=\sqrt{-{\displaystyle \frac{{\nu }_2}{{\nu }_4}}}sech\left(\sqrt{{\nu }_2}{\xi }_2\right)$$

  (45)

  for ${\nu }_2>0,{\nu }_4<0$. The wave coordinates ${\xi }_1$ and ${\xi }_2$ are given in Equation (37).
• Family 2: A solution combined distinct generalized trigonometric functions with different independent variables:

  $$\begin{array}{ccc}& & u_1^{\left(2\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(2\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (46)

  where

  $$V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(2\right)}({\xi }_1;1,2,4)=\sqrt{{\displaystyle \frac{{\mu }_2}{2{\mu }_4}}}tan\left({\displaystyle \frac{{\mu }_2}{2}}{\xi }_1\right)$$

  (47)

  for ${\mu }_0={\displaystyle \frac{{\mu }_2^2}{4{\mu }_4}}$, ${\mu }_2>0,{\mu }_4<0$ and

  $$V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(2\right)}({\xi }_2;1,2,4)=\sqrt{-{\displaystyle \frac{{\nu }_2}{{\nu }_4}}}sec\left(\sqrt{-{\nu }_2}{\xi }_2\right)$$

  (48)

  for ${\nu }_2<0,{\nu }_4>0$. The wave coordinates ${\xi }_1$ and ${\xi }_2$ are given in Equation (37).
• Family 3: Solutions combined generalized hyperbolic and trigonometric functions with different independent variables:

  $$\begin{array}{ccc}& & u_1^{\left(3\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(3\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (49)

  where ${V^{\left(1\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}({\xi }_1;1,2,4)$ is presented in Equation (44) and $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(2\right)}({\xi }_1;1,2,4)$ is presented in Equation (48). The wave coordinates ${\xi }_1$ and ${\xi }_2$ are given in Equation (37).

  $$\begin{array}{ccc}& & u_1^{\left(4\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(4\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (50)

  where ${V^{\left(2\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}({\xi }_1;1,2,4)$ is presented in Equation (47) and $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(2\right)}({\xi }_1;1,2,4)$ is presented in Equation (45). The wave coordinates ${\xi }_1$ and ${\xi }_2$ are given in Equation (37).
• Family 4: Solutions combined generalized elliptic and hyperbolic functions with different independent variables:

  $$\begin{array}{ccc}& & u_1^{\left(5\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(3\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(5\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(3\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (51)

  where

  $$V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(3\right)}({\xi }_1;1,2,4)=\sqrt{-{\displaystyle \frac{{\mu }_2{m}^2}{{\mu }_4(2m^2-1)}}}cn\left(\sqrt{{\displaystyle \frac{{\mu }_2{m}^2}{2m^2-1}}}{\xi }_1\right)$$

  (52)

  for ${\mu }_0={\displaystyle \frac{{\mu }_2^2{m}^2(m^2-1)}{{\mu }_4{(2m^2-1)}^2}},{\mu }_2>0,{\mu }_4<0$. The special function $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(1\right)}({\xi }_2;1,2,4)$ is presented in Equation (45).

  $$\begin{array}{ccc}& & u_1^{\left(6\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(4\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(6\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(4\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (53)

  where

  $$V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(4\right)}({\xi }_1;1,2,4)=\sqrt{-{\displaystyle \frac{{\mu }_2{m}^2}{{\mu }_4(2-m^2)}}}dn\left(\sqrt{{\displaystyle \frac{{\mu }_2{m}^2}{2-m^2}}}{\xi }_1\right)$$

  (54)

  for ${\mu }_0={\displaystyle \frac{{\mu }_2^2(1-m^2)}{{\mu }_4{(m^2-2)}^2}},{\mu }_2>0,{\mu }_4<0$.The special function $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(1\right)}({\xi }_2;1,2,4)$ is presented in Equation (45).

  $$\begin{array}{ccc}& & u_1^{\left(7\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(5\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(7\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(5\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(1\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (55)

  where

  $$V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(5\right)}({\xi }_1;1,2,4)=\pm \sqrt{-{\displaystyle \frac{{\mu }_2{m}^2}{{\mu }_4(m^2+1)}}}sn\left(\sqrt{{\displaystyle \frac{-{\mu }_2}{m^2+1}}}{\xi }_1\right)$$

  (56)

  for ${\mu }_0={\displaystyle \frac{{\mu }_2^2{m}^2}{{\mu }_4{(m^2+1)}^2}},{\mu }_2<0,{\mu }_4>0$.The special function $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(1\right)}({\xi }_2;1,2,4)$ is presented in Equation (45). In all solutions in this paragraph, the wave coordinates ${\xi }_1$ and ${\xi }_2$ are presented by Equation (37).
• Family 5: Solutions combined generalized elliptic and trigonometric functions with different independent variables:

  $$\begin{array}{ccc}& & u_1^{\left(8\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(3\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(8\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(3\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (57)

  where $V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(3\right)}({\xi }_1;1,2,4)$ is presented in Equation (52) and $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(2\right)}({\xi }_2;1,2,4)$ is presented in Equation (48).

  $$\begin{array}{ccc}& & u_1^{\left(9\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(4\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(9\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(4\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (58)

  where $V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(4\right)}({\xi }_1;1,2,4)$ is presented in Equation (54) and $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(2\right)}({\xi }_2;1,2,4)$ is presented in Equation (48).

  $$\begin{array}{ccc}& & u_1^{\left(10\right)}({\xi }_1,{\xi }_2)={ϵ}_0+{\displaystyle \frac{3{\mu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(5\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\eta }_0-{\displaystyle \frac{3{\nu }_4}{2}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & v_1^{\left(10\right)}({\xi }_1,{\xi }_2)=-{\displaystyle \frac{3{\mu }_4}{4}}\left({\displaystyle \frac{2-2{\eta }_0-{\omega }_2-{\omega }_1-2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2}{{\mu }_2-{\nu }_2}}\right){V^{\left(5\right)}}_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \\ & & +{\displaystyle \frac{1}{2}}{\eta }_0+{\displaystyle \frac{1}{4}}{\omega }_1-5\alpha {\mu }_2+{\displaystyle \frac{1}{4}}{\omega }_2-5\alpha {\nu }_2-{\theta }_0+{\displaystyle \frac{1}{2}}{ϵ}_0\hfill \\ & & -{\displaystyle \frac{3{\nu }_4}{4}}\left({\displaystyle \frac{2{\eta }_0+{\omega }_2+{\omega }_1+2{ϵ}_0+20\alpha {\mu }_2-20\alpha {\nu }_2-2}{{\mu }_2-{\nu }_2}}\right){V^{\left(2\right)}}_{0,0,{\mu }_2,0,{\mu }_4}^2({\xi }_1;1,2,4)\hfill \end{array}$$

  (59)

  where $V_{{\mu }_0,0,{\mu }_2,0,{\mu }_4}^{\left(5\right)}({\xi }_1;1,2,4)$ is presented in Equation (56) and $V_{0,0,{\mu }_2,0,{\mu }_4}^{\left(2\right)}({\xi }_2;1,2,4)$ is presented in Equation (48). In all the solutions in this paragraph, the wave coordinates ${\xi }_1$ and ${\xi }_2$ are are presented by Equation (37).

Illustrative examples of two exact solutions of Equation (1) are given in f1fig,f2fig.

Remark 1: The hyperbolic, trigonometric and elliptic functions that appear in the solutions of this section are defined as ’generalized` functions with respect to their arguments. It is evident that these arguments include fractional order terms, such as ${\displaystyle \frac{t^{\varphi }}{\Gamma (1+\varphi )}}$, which reflect the memory and nonlocal effects inherent in fractional calculus. It is obvious that the generalized functions thus defined can easily be reduced to the standard functions when the fractional order $\varphi =1$, thus restoring the classical (integer–order) case. Therefore, for the particular case $\varphi =1$, all exact solutions presented above also satisfy the classical (integer–order) version of the system (1).

Remark 2: When ${\omega }_1={\omega }_2$, the exact solutions presented above reduce to those discussed in Step 1 (variant 2b) of the generalized SEsM algorithm, i.e. the distinct types of special functions present in the obtained analytical solutions have the same independent variable. To reach this particular case, it is easy to substitute ${\omega }_1={\omega }_2=\omega$ in Equation (42) and to rewrite all special functions present in the exact solutions of Equation (1) at ${\xi }_1={\xi }_2=\xi$.

Figure 1. Numerical simulations of the solution (50) for ${ϵ}_0=20,{\eta }_0=0.4,{\mu }_2=0.5,{\mu }_4=-2;{\nu }_2=0.2,{\nu }_4=-0.5;\alpha =2,{\omega }_1=0.05,{\omega }_2=0.1$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

Figure 1. Numerical simulations of the solution (50) for ${ϵ}_0=20,{\eta }_0=0.4,{\mu }_2=0.5,{\mu }_4=-2;{\nu }_2=0.2,{\nu }_4=-0.5;\alpha =2,{\omega }_1=0.05,{\omega }_2=0.1$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

Figure 2. Numerical simulations of the solution (55) for ${ϵ}_0=20,{\eta }_0=0.4,{\mu }_2=-5,{\mu }_4=2;{\nu }_2=0.2,{\nu }_4=-0.4;\alpha =2,{\omega }_1=0.5,{\omega }_2=0.1$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

Figure 2. Numerical simulations of the solution (55) for ${ϵ}_0=20,{\eta }_0=0.4,{\mu }_2=-5,{\mu }_4=2;{\nu }_2=0.2,{\nu }_4=-0.4;\alpha =2,{\omega }_1=0.5,{\omega }_2=0.1$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

5. Derivation of Multi–Wave Exact Solutions of Equation (2) Using a Standard Traveling Wave Transformation

We present the general solution of Equation (2) again as combination of two single composite functions including two independent variables:

$$\begin{array}{ccc}& & u_2({\xi }_1,{\xi }_2)=\sum _{i_1=0}^{n_1}{ϵ}_{i_1}{\left[f_1\left({\xi }_1\right)\right]}^{i_1}+\sum _{i_2=0}^{n_2}{\eta }_{i_2}{\left[f_2\left({\xi }_2\right)\right]}^{i_2}\hfill \\ & & v_2({\xi }_1,{\xi }_2)=\sum _{i_1=0}^{n_3}{\theta }_{i_1}{\left[f_1\left({\xi }_1\right)\right]}^{i_1}+\sum _{i_2=0}^{n_4}{\iota }_{i_2}{\left[f_2\left({\xi }_2\right)\right]}^{i_2}\hfill \end{array}$$

(60)

where the following standart traveling wave transformations are introfuced:

$${\xi }_1=x-{\omega }_{1}t,{\xi }_2=x-{\omega }_{2}t$$

(61)

We shall use the fractional variants of the ordinary differential equations of Riccati as simple equations:

$$D_{{\xi }_1}^{\varphi }{f}_1={\mu }_0+f_1^2,D_{{\xi }_2}^{\varphi }{f}_2={\nu }_0+f_2^2,0\le \varphi \le 1,$$

(62)

The analytical solutions of Equation (62) are given in Appendix A.2 and are defined in sense of the modified Jumarie’s Riemann-Liouville derivatives.

The balance equations derived are $n_1=n_3=2m_1-2$ and $n_2=n_4=2m_2-2$. Bellow we shall present the solution (60) again by special functions $V_{{\mu }_0,0,1}({\xi }_1,\varphi ;2)$ and $V_{{\nu }_0,0,1}({\xi }_2,\varphi ;2)$, as follows:

$$\begin{array}{ccc}& & u_2({\xi }_1,{\xi }_2)={ϵ}_0+{ϵ}_1{V}_{{\mu }_0,0,1}({\xi }_1,\varphi ;2)+{ϵ}_2{V}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\eta }_0+{\eta }_1{V}_{{\nu }_0,0,1}({\xi }_2,\varphi ;2)+{\eta }_2{V}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2),\hfill \\ & & v_2({\xi }_1,{\xi }_2)={\theta }_0+{\theta }_1{V}_{{\mu }_0,0,1}({\xi }_1,\varphi ;2)+{\theta }_2{V}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+{\iota }_1{V}_{{\nu }_0,0,1}({\xi }_2,\varphi ;2)++{\iota }_2{V}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

(63)

Substitution of Equation (60), Equation (61) and Equation (62) in Equation (2) leads to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}& & 8{ϵ}_2{\theta }_1+8{ϵ}_2{\theta }_2-48k{ϵ}_1-48k{ϵ}_2+8{ϵ}_1{ϵ}_2+8{ϵ}_1{\theta }_1+8{ϵ}_1{\theta }_2+4{{ϵ}_1}^2+4{{ϵ}_2}^2=0\hfill \\ & & 4{\eta }_2{ϵ}_1+4{\eta }_2{ϵ}_2+4{\eta }_2{\theta }_1+4{\eta }_2{\theta }_2+4{\iota }_2{ϵ}_1+4{\iota }_2{ϵ}_2=0\hfill \\ & & 4{\eta }_2{ϵ}_1+4{\eta }_2{ϵ}_2+4{\eta }_2{\theta }_1+4{\eta }_2{\theta }_2+4{\iota }_2{ϵ}_1+4{\iota }_2{ϵ}_2=0\hfill \\ & & 8{ϵ}_2{\mu }_0{\theta }_2+8{ϵ}_1{\mu }_0{\theta }_1+8{ϵ}_1{\mu }_0{ϵ}_2+8{ϵ}_2{\mu }_0{\theta }_1\hfill \\ & & +8{ϵ}_1{\mu }_0{\theta }_2+2{\omega }_1^{\varphi }{\theta }_1-2{\omega }_1^{\varphi }{ϵ}_1-2{\omega }_1{ϵ}_2-2{\omega }_2^{\varphi }{ϵ}_1-2{\omega }_2^{\varphi }{ϵ}_2+4{ϵ}_0{ϵ}_1\hfill \\ & & +4{ϵ}_0{ϵ}_2+4{{ϵ}_1}^2{\mu }_0+4{{ϵ}_2}^2{\mu }_0+4{\eta }_0{ϵ}_1+4{\eta }_0{ϵ}_2\hfill \\ & & +4{\theta }_0{ϵ}_1+4{\theta }_0{ϵ}_2+4{\iota }_0{ϵ}_1+4{\iota }_0{ϵ}_2+2{\omega }_1^{\varphi }{\theta }_2+2{\omega }_2^{\varphi }{\theta }_1+2{\omega }_2^{\varphi }{\theta }_2\hfill \\ & & +4{ϵ}_0{\theta }_1+4{ϵ}_0{\theta }_2+4{\eta }_0{\theta }_1+4{\eta }_0{\theta }_2=0\hfill \\ & & 2{\eta }_1{ϵ}_1+2{\iota }_1{ϵ}_1+2{\eta }_1{ϵ}_2+2{\eta }_1{\theta }_2+2{\iota }_1{ϵ}_2+2{\eta }_1{\theta }_1=0\hfill \\ & & 4{\theta }_2{\eta }_2{\nu }_0+4{ϵ}_2{\iota }_2{\nu }_0+4{\theta }_1{\eta }_2{\nu }_0+4{ϵ}_2{\eta }_2{\nu }_0+4{ϵ}_1{\eta }_2{\nu }_0+4{ϵ}_1{\iota }_2{\nu }_0=0\hfill \\ & & -2{\omega }_2^{\varphi }{\mu }_0{ϵ}_2+4{ϵ}_0{\mu }_0{ϵ}_1+4{ϵ}_0{\mu }_0{ϵ}_2+4{\eta }_0{\mu }_0{ϵ}_1\hfill \\ & & +4{\eta }_0{\mu }_0{ϵ}_2+2{\omega }_1{\mu }_0{\theta }_1+2{\omega }_1^{\varphi }{\mu }_0{\theta }_2+2{\omega }_2^{\varphi }{\mu }_0{\theta }_1\hfill \\ & & +2{\omega }_2^{\varphi }{\mu }_0{\theta }_2+4{ϵ}_0{\mu }_0{\theta }_1+4{ϵ}_0{\mu }_0{\theta }_2-2{\omega }_1^{\varphi }{\mu }_0{ϵ}_1-2{\omega }_1^{\varphi }{\mu }_0{ϵ}_2\hfill \\ & & -2{\omega }_2^{\varphi }{\mu }_0{ϵ}_1+4{\mu }_0{\theta }_1+4{\mu }_0{\theta }_2+4{\eta }_0{\mu }_0{\theta }_1+4{\eta }_0{\mu }_0{\theta }_2\hfill \\ & & +4{\theta }_0{\mu }_0{ϵ}_1+4{\theta }_0{\mu }_0{ϵ}_2+4{\iota }_0{\mu }_0{ϵ}_1+4{\iota }_0{\mu }_0{ϵ}_2-32{{\mu }_0}^2{ϵ}_1-32{{\mu }_0}^2{ϵ}_2=0\hfill \\ & & 8{\eta }_2{\iota }_2-48{\eta }_2+4{{\eta }_2}^2=0\hfill \\ & & 6{\eta }_1{\eta }_2+6{\eta }_2{\iota }_1-12{\eta }_1+6{\eta }_1{\iota }_2=0\hfill \\ & & 8{\eta }_2{\iota }_2{\nu }_0+4{\iota }_2+4{ϵ}_0{\iota }_2+4{{\eta }_2}^2{\nu }_0+2{\omega }_1^{\varphi }{\iota }_2-2{\omega }_2{\eta }_2-80{\eta }_2{\nu }_0\hfill \\ & & +2{\omega }_2^{\varphi }{\iota }_2+2{{\eta }_1}^2+4{\eta }_0{\iota }_2+4{\iota }_0{\eta }_2+4{\theta }_0{\eta }_2+4{\eta }_1{\iota }_1+4{\eta }_0{\eta }_2+4{ϵ}_0{\eta }_2-2{\omega }_1^{\varphi }{\eta }_2=0\hfill \\ & & {\omega }_1^{\varphi }{\iota }_1{\nu }_0-{\omega }_2{\eta }_1{\nu }_0+{\omega }_2{\iota }_1{\nu }_0+2{\eta }_0{\eta }_1{\nu }_0+2{ϵ}_0{\eta }_1{\nu }_0\hfill \\ & & +2{ϵ}_0{\iota }_1{\nu }_0+2{\iota }_1{\nu }_0-{\omega }_1^{\varphi }{\eta }_1{\nu }_0+2{\eta }_0{\iota }_1{\nu }_0-4{\eta }_1{{\nu }_0}^2+2{\iota }_0{\eta }_1{\nu }_0+2{\theta }_0{\eta }_1{\nu }_0=0\hfill \end{array}$$

(64)

There are many non–trivial solutions of this algebraic system. For example, one simpler non–trivial solution is:

$$\begin{array}{ccc}& & {ϵ}_1={\eta }_1={\theta }_1={\iota }_1=0,{ϵ}_2={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}},\hfill \\ & & {\eta }_0=8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_1^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1,{\eta }_2=-{\displaystyle \frac{1}{2}}{ϵ}_2\hfill \\ & & {\theta }_2=-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}},{\iota }_2=-2{\theta }_2,{ϵ}_0=0\hfill \end{array}$$

(65)

where ${\theta }_0,{\iota }_0,{\omega }_1,{\omega }_2,{\mu }_0,{\nu }_0$ are free parameters. Substituting the coefficients of Equation (65) in Equation (63) and combining the different analytical solution of the second equation in (39) leads to the following families of exact solutions of Equation (2):

• Family 1: Solutions combined distinct fractional generalized hyperbolic functions with different independent variables (When ${\mu }_0<0$ and ${\nu }_0<0$):

  $$\begin{array}{ccc}& & u_2^{\left(1\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_1^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(1\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (66)

  where

  $$V_{{\mu }_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)=-\sqrt{-{\mu }_0}tanh[\sqrt{-{\mu }_0}{\xi }_1,\varphi ]$$

  (67)

  $$V_{{\nu }_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)=-\sqrt{-{\nu }_0}tanh[\sqrt{-{\nu }_0}{\xi }_2,\varphi ]$$

  (68)

  $$\begin{array}{ccc}& & u_2^{\left(2\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_1^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(2\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (69)

  where

  $$V_{{\mu }_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2)=-\sqrt{-{\mu }_0}coth[\sqrt{-{\mu }_0}{\xi }_1,\varphi ]$$

  (70)

  $$V_{{\nu }_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)=-\sqrt{-{\nu }_0}coth[\sqrt{-{\nu }_0}{\xi }_2,\varphi ]$$

  (71)

  $$\begin{array}{ccc}& & u_2^{\left(3\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_1^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(3\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (72)

  where$V_{{\mu }_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (67) and $V_{{\nu }_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (71).

  $$\begin{array}{ccc}& & u_2^{\left(4\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_1^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(4\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (73)

  where$V_{{\mu }_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;1,2)$ is presented in Equation (70) and $V_{{\nu }_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (68).
• Family 2: Solutions combined distinct fractional generalized trigonometric functions with different independent variables (When ${\mu }_0>0$ and ${\nu }_0>0$):

  $$\begin{array}{ccc}& & u_2^{\left(5\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(1\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (74)

  where

  $$V_{{\mu }_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)=\sqrt{{\mu }_0}tan[\sqrt{{\mu }_0}{\xi }_1,\varphi ]$$

  (75)

  $$V_{{\nu }_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)=\sqrt{{\nu }_0}tan[\sqrt{{\nu }_0}{\xi }_2,\varphi ]$$

  (76)

  $$\begin{array}{ccc}& & u_2^{\left(6\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(6\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (77)

  where

  $$V_{{\mu }_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)=-\sqrt{-{\mu }_0}cot[\sqrt{{\mu }_0}{\xi }_1,\varphi ]$$

  (78)

  $$V_{{\nu }_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)=-\sqrt{-{\nu }_0}cot[\sqrt{{\nu }_0}{\xi }_2,\varphi ]$$

  (79)

  $$\begin{array}{ccc}& & u_2^{\left(7\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(7\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (80)

  where $V_{{\mu }_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (75) and $V_{{\nu }_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (79).

  $$\begin{array}{ccc}& & u_2^{\left(8\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(8\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (81)

  where $V_{{\mu }_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (78 ) and $V_{{\nu }_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (76).
• Family 3: Solutions combined distinct fractional generalized hyperbolic and trigonometric functions with different independent variables
  • (a) When ${\mu }_0<0$ and ${\nu }_0>0$

    $$\begin{array}{ccc}& & u_2^{\left(9\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(9\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (82)

    where $V_{{\mu }_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (67) and $V_{{\nu }_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (76).

    $$\begin{array}{ccc}& & u_2^{\left(10\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(10\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (83)

    where $V_{{\mu }_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (67) and $V_{{\nu }_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (79).

    $$\begin{array}{ccc}& & u_2^{\left(11\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(11\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (84)

    where $V_{{\mu }_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (70) and $V_{{\nu }_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (76).

    $$\begin{array}{ccc}& & u_2^{\left(12\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(12\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (85)

    where $V_{{\mu }_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (70) and $V_{{\nu }_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (79).
  • (b) When ${\mu }_0>0$ and ${\nu }_0<0$

    $$\begin{array}{ccc}& & u_2^{\left(13\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(13\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (86)

    where $V_{{\mu }_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (75) and $V_{{\nu }_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (68).

    $$\begin{array}{ccc}& & u_2^{\left(14\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(14\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(3\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (87)

    where $V_{{\mu }_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (75) and $V_{{\nu }_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (71).

    $$\begin{array}{ccc}& & u_2^{\left(15\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(15\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(1\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (88)

    where $V_{{\mu }_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)$ is presented in in Equation (78) and $V_{{\nu }_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)$ is presented in in Equation (67).

    $$\begin{array}{ccc}& & u_2^{\left(16\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(16\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(4\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(2\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

    (89)

    where $V_{{\mu }_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)$ is presented in Equation (78) and $V_{{\nu }_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)$ is presented in Equation (71).
• Familly 3. Solutions combained distinct fractional generalized rational (algebraic) functions with different independent variables (When ${\mu }_{00}=0$ and ${\nu }_0=0$)

  $$\begin{array}{ccc}& & u_2^{\left(17\right)}({\xi }_1,{\xi }_2)={\displaystyle \frac{3(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(5\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +8{\mu }_0+8{\nu }_0-{ϵ}_0-2{\theta }_0-2{\iota }_0+{\displaystyle \frac{3}{2}}{\omega }_5^{\varphi }+{\displaystyle \frac{3}{2}}{\omega }_2^{\varphi }+1\hfill \\ & & -{\displaystyle \frac{3}{2}}{\displaystyle \frac{(1+{\omega }_2^{\varphi }+4{\mu }_0+4{\nu }_0-{\theta }_0-{\iota }_0+{\omega }_1^{\varphi })}{{\mu }_0-{\nu }_0}}{V^{\left(5\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \\ & & v_2^{\left(15\right)}({\xi }_1,{\xi }_2)={\theta }_0-{\displaystyle \frac{3}{2}}{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(5\right)}}_{{\mu }_0,0,1}^2({\xi }_1,\varphi ;2)\hfill \\ & & +{\iota }_0+3{\displaystyle \frac{1-{\theta }_0+{\omega }_2^{\varphi }+{\omega }_1^{\varphi }-{\iota }_0+8{\nu }_0}{{\mu }_0-{\nu }_0}}{V^{\left(5\right)}}_{{\nu }_0,0,1}^2({\xi }_2,\varphi ;2)\hfill \end{array}$$

  (90)

  where

  $$V_{{\mu }_0,0,1}^{\left(5\right)}({\xi }_1,\varphi ;2)={\displaystyle \frac{\Gamma (1+\varphi )}{{\xi }_1^{\varphi }+{\omega }_1}}$$

  (91)

  $$V_{{\nu }_0,0,1}^{\left(5\right)}({\xi }_2,\varphi ;2)={\displaystyle \frac{\Gamma (1+\varphi )}{{\xi }_2^{\varphi }+{\omega }_2}}$$

  (92)

  In all the solutions provided in this section, the wave coordinates ${\xi }_1$ and${\xi }_2$ are presented by Equation (61)

The fractional generalized trigonometric and hyperbolic functions presented in the most solutions are defined as follows:

$$\begin{array}{ccc}& & sin({\xi }_1,\varphi )={\displaystyle \frac{E_{\varphi }\left(i{\xi }_1^{\varphi }\right)-E_{\varphi }(-i{\xi }_1^{\varphi })}{2i}},sin({\xi }_2,\varphi )={\displaystyle \frac{E_{\varphi }\left(i{\xi }_2^{\varphi }\right)-E_{\varphi }(-i{\xi }_2^{\varphi })}{2i}},\hfill \\ & & cos({\xi }_1,\varphi )={\displaystyle \frac{E_{\varphi }\left(i{\xi }_1^{\varphi }\right)+E_{\varphi }(-i{\xi }_1^{\varphi })}{2i}},cos({\xi }_2,\varphi )={\displaystyle \frac{E_{\varphi }\left(i{\xi }_2^{\varphi }\right)+E_{\varphi }(-i{\xi }_2^{\varphi })}{2i}},\hfill \\ & & tan({\xi }_1,\varphi )={\displaystyle \frac{sin({\xi }_1,\varphi )}{cos({\xi }_1,\varphi )}},tan({\xi }_2,\varphi )={\displaystyle \frac{sin({\xi }_2,\varphi )}{cos({\xi }_2,\varphi )}},\hfill \\ & & cot({\xi }_1,\varphi )={\displaystyle \frac{cos({\xi }_1,\varphi )}{sin({\xi }_1,\varphi )}},cot({\xi }_2,\varphi )={\displaystyle \frac{cos({\xi }_2,\varphi )}{sin({\xi }_2,\varphi )}},\hfill \\ & & sinh({\xi }_1,\varphi )={\displaystyle \frac{E_{\varphi }\left({\xi }_1^{\varphi }\right)-E_{\varphi }\left({\xi }_1^{\varphi }\right)}{2}},sinh({\xi }_2,\varphi )={\displaystyle \frac{E_{\varphi }\left({\xi }_2^{\varphi }\right)-E_{\varphi }\left({\xi }_2^{\varphi }\right)}{2}},\hfill \\ & & cosh({\xi }_1,\varphi )={\displaystyle \frac{E_{\varphi }\left({\xi }_1^{\varphi }\right)+E_{\varphi }\left({\xi }_1^{\varphi }\right)}{2}},cosh({\xi }_2,\varphi )={\displaystyle \frac{E_{\varphi }\left({\xi }_2^{\varphi }\right)+E_{\varphi }\left({\xi }_2^{\varphi }\right)}{2}},\hfill \\ & & tanh({\xi }_1,\varphi )={\displaystyle \frac{sinh({\xi }_1,\varphi )}{cosh({\xi }_1,\varphi )}},tanh({\xi }_2,\varphi )={\displaystyle \frac{sinh({\xi }_2,\varphi )}{cosh({\xi }_2,\varphi )}},\hfill \\ & & coth({\xi }_1,\varphi )={\displaystyle \frac{cosh({\xi }_1,\varphi )}{sinh({\xi }_1,\varphi )}},coth({\xi }_2,\varphi )={\displaystyle \frac{cosh({\xi }_2,\varphi )}{sinh({\xi }_2,\varphi )}},\hfill \end{array}$$

(93)

where $E_{\varphi }\left({\xi }_1\right)={\sum }_{k=0}^{\infty }{\displaystyle \frac{{\xi }_1^k}{\Gamma (1+k\varphi )}}$ and $E_{\varphi }\left({\xi }_2\right)={\sum }_{k=0}^{\infty }{\displaystyle \frac{{\xi }_2^k}{\Gamma (1+k\varphi )}}$ are Mittag–Leffler functions.

Illustrative examples of two exact solutions of Equation (1) are given in f3fig,f4fig.

Remark 1: The hyperbolic, trigonometric, and rational functions that appear in the solutions of this section are defined as “fractional generalized” functions with respect to the fractional parameter $\varphi \in (0,1]$. They differ from the generalized hyperbolic and trigonometric functions of the form that were introduced in the previous section, but reduce back to their corresponding classical counterparts at $\varphi =1$. It is obvious that the fractional generalized functions so defined can easily be reduced to the standard functions when the fractional order $\varphi =1$, thus restoring the classical (integer) case. Therefore, for the particular case $\varphi =1$, all the exact solutions presented above also are valid the classical (integer) version of the system (2) at $\varphi =1$.

Remark 2: When ${\omega }_1={\omega }_2$, the exact solutions presented above reduce to those discussed in Step 1 (variant 2b) of the generalized SEsM algorithm, i.e. the distinct types of special functions present in the obtained analytical solutions have the same independent variable. To achieve this particular case, it is easy to replace ${\omega }_1={\omega }_2=\omega$ in Equation (65) and rewrite all special functions present in the exact solutions of Equation (1) at ${\xi }_1={\xi }_2=\xi$.

Figure 3. Numerical simulations of the solution (80) for ${\iota }_0=20,{\theta }_0=0.4,{\mu }_0=50,{\nu }_0=0.2,{\omega }_1=0.5,{\omega }_2=0.1$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

Figure 3. Numerical simulations of the solution (80) for ${\iota }_0=20,{\theta }_0=0.4,{\mu }_0=50,{\nu }_0=0.2,{\omega }_1=0.5,{\omega }_2=0.1$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

Figure 4. Numerical simulations of the solution (83) for ${\iota }_0=20,{\theta }_0=0.4,{\mu }_0=50,{\nu }_0=0.2,{\omega }_1=0.05,{\omega }_2=0.01$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

Figure 4. Numerical simulations of the solution (83) for ${\iota }_0=20,{\theta }_0=0.4,{\mu }_0=50,{\nu }_0=0.2,{\omega }_1=0.05,{\omega }_2=0.01$: 3D plots of (a) u(x,t) and (b) v(x,t); 2D plots of (c) u(x,t) and (d) v(x,t) for $t=1$

6. Conclusions

In this paper, a generalized algorithm based on the extended SEsM is proposed for obtaining exact solutions of FNPDEs (in particular, of NPDEs). This algorithm provides several options at each step, which allows the user to choose the most appropriate option depending on the expected wave dynamics in the modeled physical context. One of the main innovations in the proposed algorithm is the ability to construct the exact solutions of the studied FNPDEs by means of composite functions, including solutions of at least two simple equations with different independent variables. Multiple configurations for constructing analytical solutions of systems of type (19) are also shown. The algorithm also includes options for choosing an appropriate wave transform, as well as an appropriate type of ODEs in a role of simple equations. In this sense, the SEsM algorithm proposed appears as a generalization (with an extension) of all ansatz–based method algorithms for finding exact solutions of generalized models described by FNPDEs (and as a particular case, by NPDEs).

To demonstrate the wide possibilities provided by the generalized SEsM algorithm, two examples of its application to two variants of a system of NPDEs, which are generalized models of wave dynamics in shallow water processes, are given. To generate exact solutions to both considered systems, it is assumed that each system variable in the studied models exhibits multi–wave behavior, which is expressed as a superposition of two waves propagating at different velocities. In extracting the exact solutions to each studied system of FNPDEs, different wave transformations and different types of simple equations were used. As a result, many new multi–wave solutions were derived, including combinations of generalized hyperbolic, trigonometric, elliptic and rational functions. Numerical illustrations of some of the obtained analytical solutions confirm the analytically presented multi–wave dynamics, expressed in superposition and overlay of various multi–soliton-like and periodic-like waves and packets. It is shown that with minimal mathematical tricks, the analytical solutions obtained in this paper can describe the complex wave dynamics of classical (integer) shallow water models. In this sense, the obtained analytical solutions have a potential application to various scientific fields such as fluid dynamics, geophysics, oceanography and others.

Appendix A.1. Analytical Solutions of the Simple Equations Used in sec4

The equation ${\left({\displaystyle \frac{df}{d\xi }}\right)}^2={\mu }_0+{\mu }_2{f}^2+{\mu }_4{f}^4$ admits the following solutions [51]:

• a hyperbolic solution
  $f\left(\xi \right)=\sqrt{{\displaystyle \frac{-{\mu }_2}{2{\mu }_4}}}tanh\left({\displaystyle \frac{{\mu }_2}{2}}\xi \right)$ for ${\mu }_0={\displaystyle \frac{{\mu }_2^2}{4{\mu }_4}}$ and ${\mu }_2<0,{\mu }_4>0$;
• a trigonometric solution
  $f\left(\xi \right)=\sqrt{{\displaystyle \frac{{\mu }_2}{2{\mu }_4}}}tan\left({\displaystyle \frac{{\mu }_2}{2}}\xi \right)$ for ${\mu }_0={\displaystyle \frac{{\mu }_2^2}{4{\mu }_4}}$ and ${\mu }_2>0,{\mu }_4<0$;
• Jacobi elliptic function solutions:
  $f\left(\xi \right)=\sqrt{-{\displaystyle \frac{{\mu }_2{m}^2}{{\mu }_4(2m^2-1)}}}cn\left(\sqrt{{\displaystyle \frac{{\mu }_2{m}^2}{2m^2-1}}}\xi \right)$ for ${\mu }_0={\displaystyle \frac{{\mu }_2^2{m}^2(m^2-1)}{{\mu }_4{(2m^2-1)}^2}},{\mu }_2>0,{\mu }_4<0$;
  $f\left(\xi \right)=\sqrt{-{\displaystyle \frac{{\mu }_2{m}^2}{{\mu }_4(2-m^2)}}}dn\left(\sqrt{{\displaystyle \frac{{\mu }_2{m}^2}{2-m^2}}}\xi \right)$ for ${\mu }_0={\displaystyle \frac{{\mu }_2^2(1-m^2)}{{\mu }_4{(m^2-2)}^2}},{\mu }_2>0,{\mu }_4<0$;
  $f\left(\xi \right)=\pm \sqrt{-{\displaystyle \frac{{\mu }_2{m}^2}{{\mu }_4(m^2+1)}}}sn\left(\sqrt{{\displaystyle \frac{-{\mu }_2}{m^2+1}}}\xi \right)$ for ${\mu }_0={\displaystyle \frac{{\mu }_2^2{m}^2}{{\mu }_4{(m^2+1)}^2}},{\mu }_2<0,{\mu }_4>0$
  where m is a modulus. The Jacobi elliptic functions are doubly periodical and possess properties of triangular functions. Additionally, we note that when $m\to 1$, the Jacobi functions degenerate to the hyperbolic functions, i.e.

  $$sn\left(\xi \right)\to tanh\left(\xi \right),cn\left(\xi \right)\to sech\left(\xi \right),dn\left(\xi \right)\to sech\left(\xi \right)$$

  When $m\to 0$, the Jacobi functions degenerate to the triangular functions, i.e.

  $$sn\left(\xi \right)\to sin\left(\xi \right),cn\left(\xi \right)\to cos\left(\xi \right),dn\left(\xi \right)\to 1$$

  Additional sub–variants of the elliptic solutions can be found in [58].

The equation ${\left({\displaystyle \frac{df}{d\xi }}\right)}^2={\nu }_2{f}^2+{\nu }_4{f}^4$ admits the following solutions [51]:

• a hyperbolic solution
  $f\left(\xi \right)=\sqrt{-{\displaystyle \frac{{\nu }_2}{{\nu }_4}}}sech\left(\sqrt{{\nu }_2}\xi \right)$ for ${\nu }_2>0,{\nu }_4<0$;
• a trigonometric solution
  $f\left(\xi \right)=\sqrt{{\displaystyle \frac{-{\nu }_2}{{\nu }_4}}}sec\left(\sqrt{-{\nu }_2}\xi \right)$ for ${\nu }_2<0,{\nu }_4>0$.

Appendix A.2. Analytical Solutions of the Simple Equations Used in sec5

The equation $D_{\xi }^{\varphi }={\mu }_0+f^2$ (or $D_{\xi }^{\varphi }={\nu }_0+f^2$) admits the following solutions [40]:

• $f\left(\xi \right)=-\sqrt{-{\mu }_0}{tanh}_{\varphi }\left(\sqrt{-{\mu }_0}\xi \right)$ for ${\mu }_0<0$;
• $f\left(\xi \right)=-\sqrt{-{\mu }_0}{coth}_{\varphi }\left(\sqrt{-{\mu }_0}\xi \right)$ for ${\mu }_0<0$;
• $f\left(\xi \right)=\sqrt{{\mu }_0}{tan}_{\varphi }\left(\sqrt{{\mu }_0}\xi \right)$ for ${\mu }_0>0$;
• $f\left(\xi \right)=-\sqrt{{\mu }_0}{cot}_{\varphi }\left(\sqrt{{\mu }_0}\xi \right)$ for ${\mu }_0>0$;
• $f\left(\xi \right)=-{\displaystyle \frac{\Gamma (1+\varphi )}{{\xi }^{\varphi }+w}}$ for ${\mu }_0>0$ and $w=const.$,

where

$$\begin{array}{ccc}& & sin(\xi ,\varphi )={\displaystyle \frac{E_{\varphi }\left(i{\xi }^{\varphi }\right)-E_{\varphi }(-i{\xi }^{\varphi })}{2i}},cos(\xi ,\varphi )={\displaystyle \frac{E_{\varphi }\left(i{\xi }^{\varphi }\right)+E_{\varphi }(-i{\xi }^{\varphi })}{2i}},\hfill \\ & & tan(\xi ,\varphi )={\displaystyle \frac{sin(\xi ,\varphi )}{cos(\xi ,\varphi )}},cot(\xi ,\varphi )={\displaystyle \frac{cos(\xi ,\varphi )}{sin(\xi ,\varphi )}},\hfill \\ & & sinh(\xi ,\varphi )={\displaystyle \frac{E_{\varphi }\left({\xi }^{\varphi }\right)-E_{\varphi }\left({\xi }^{\varphi }\right)}{2}},cosh(\xi ,\varphi )={\displaystyle \frac{E_{\varphi }\left({\xi }^{\varphi }\right)+E_{\varphi }\left({\xi }^{\varphi }\right)}{2}},\hfill \\ & & tanh(\xi ,\varphi )={\displaystyle \frac{sinh(\xi ,\varphi )}{cosh(\xi ,\varphi )}},coth(\xi ,\varphi )={\displaystyle \frac{cosh(\xi ,\varphi )}{sinh(\xi ,\varphi )}},\hfill \end{array}$$

where $E_{\varphi }\left(\xi \right)={\sum }_{k=0}^{\infty }{\displaystyle \frac{{\xi }^k}{\Gamma (1+k\varphi )}}$ is a Mittag–Leffler function.