On the Extended Simple Equations Method (SEsM) and Its Application for Finding Exact Solutions of the Time-Fractional Diffusive Predator–Prey System Incorporating an Allee Effect

Abstract

In this paper, I extend the Simple Equations Method (SEsM) and adapt it to obtain exact solutions of systems of fractional nonlinear partial differential equations (FNPDEs). The novelty in the extended SEsM algorithm is that, in addition to introducing more simple equations in the construction of the solutions of the studied FNPDEs, it is assumed that the selected simple equations have different independent variables (i.e., different coordinates moving with the wave). As a consequence, nonlinear waves propagating with different wave velocities will be observed. Several scenarios of the extended SEsM are applied to the time-fractional predator–prey model under the Allee effect. Based on this, new analytical solutions are derived. Numerical simulations of some of these solutions are presented, adequately capturing the expected diverse wave dynamics of predator–prey interactions.

1. Introduction

In the last two decades, the use of time- and space-fractional nonlinear partial differential equations (FNPDEs) has become increasingly popular in the scientific world to model and study various natural processes. B This approach enables the capture of the complexities inherent in many real-world systems, such as those found in biology and ecology [1,2,3], fluid mechanics [4,5,6,7], finance and economics [8,9,10], engineering [11,12,13], and other domains. These systems are typically characterized by memory effects, spatial heterogeneity, and non-local interactions, which cannot be adequately addressed by traditional integer-order models. One direction to elucidate these non-specific dynamical effects is to derive and analyze analytical solutions of the fractional models. Finding analytical solutions of FNPDEs, however, is a challenging task due to the complexity introduced by both the nonlinearity and the fractional order of the derivatives.

Recently, two basic approaches for extracting exact solutions of FNPDEs have been used: (1) Applying fractional transformation, which allows for reducing the studied NFPDEs to integer-order nonlinear ODEs, and then using well-established techniques from the theory of integer-order NPDEs to obtain exact general and particular solutions of NPDEs. For this scenario, numerous analytical methods have been developed, based on the pioneering works that utilize ansatz techniques [14,15,16,17,18,19,20], group analysis [21,22], perturbation methods [23,24], decomposition methods [25,26], and so on. I note that the SEsM used in this study has an almost universal character, because most known methods for finding analytical solutions to nonlinear PDEs by using integer-order ODEs or some special functions depend on particular cases, as proven in [27,28,29,30]. (2) Applying standard traveling wave transformation, which allows for reducing the studied NFPDEs to simpler fractional differential equations (in most cases, to fractional ODEs) and using their known solutions in constructing the solutions of the basic NPDEs [31]. The Mittag–Leffler function usually appears in these solutions due to its suitability in describing fractional-order systems. In this category, the most popular methods are the fractional sub-equation method [32,33,34,35], which uses known solutions of the fractional ODE of Riccati and the fractional Jacobi elliptic equation method [36,37], which uses the Jacobi elliptic functions in deriving analytical solutions of the studied FPDEs.

In general, the SEsM used in this article gives an opportunity to construct the solution of the studied single NPDE (or the studied system of NPDEs) as a composite function of one or more single functions, which include power series of solutions of one or more simple equations or of one or more special functions with one independent variable [38]. The new and different emphasis to the SEsM algorithm presented in [39,40,41,42,43,44,45,46,47,48,49,50] is the assumption that the simple equations used have different variables. This new approach was first applied to find analytical solutions of the extended integer-order KDV equation in [39]. Later, in [40], the authors independently proposed a similar methodology based on the extended Kudryashov method, as it was applied to find exact solutions of an integer-order Boussinesq-like system and the integer-order shallow-water wave equation. In terms of finding exact solutions of systems of FNPDEs, a version of the extended SEsM was applied in [49,50] to dynamical models describing natural processes in fluid mechanics and ecology.

Bellow, I shall apply another version of the extended SEsM that is applicable to specific model systems of FNPDEs, including both approaches mentioned above. In detail, as will be shown in the next section of this paper, firstly, I extend the SEsM algorithm given in [51,52,53] by presenting solutions of the studied systems by distinct single composite functions with different independent variables (different wave coordinates) for each system variable. Such a presentation of the solutions gives a more realistic picture of the wave dynamics of real-world processes whose system components exhibit different wave behaviors, including that the system variables propagate with different wave velocities. Next, I adapt the extended SEsM to finding exact solutions to studied systems of FNPDEs by introducing appropriate traveling wave transformations.

In this paper, I shall give an example of application of this variant of the extended SEsM for finding analytical solutions of the dimensionless diffusive predator–prey model, which is given in [54]:

$$\begin{array}{c}\hfill D_t^{\varphi }u=u_{xx}-\beta u+(\beta +1)u^2-u^3-(1+\alpha v)uv,\\ \hfill \kappa D_t^{\varphi }v=v_{xx}+\sigma (1+\alpha v)uv-v({\mu }_1+{\mu }_{2}v+{\mu }_3{v}^2)\end{array}$$

(1)

where $u(x,t)$ and $v(x,t)$ are dimensionless sizes or densities of the prey and predator populations, respectively, $\varphi$ is the time-fractional derivative order, and $\alpha ,\beta ,\sigma ,{\mu }_1,{\mu }_2$, and ${\mu }_3$ are dimensionless parameters. System (1) describes the complex dynamical interactions between predator and prey populations that are both subject to the Allee effect. The Allee effect refers to a phenomenon in population ecology where individuals in a population have reduced fitness or reproductive success at low population densities, leading to difficulties in survival or growth. For construction of the aforementioned model, in [54], an extension and a combination of two well-known diffusive predator–prey models with integer-order derivatives, presented in [55,56], were carried out. Subsequently, fractional derivatives in both time and space were introduced into the extended model. A traveling wave solution of the time- and space-fractional predator–prey system was derived in [54] using the extended SEsM, and the ODEs of Riccati were used as simple equations.

Here, I present briefly the basic assumptions made in constructing the time-fractional predator–prey system given in (1):

–

Like [55], in the reaction part of the first equation of (1), the local growth of the prey population incorporating the Allee effect is included, where $\beta$ is the Allee threshold below which the prey population may decline to extinction.

–

As in [56], the functional response of the predator is assumed to depend on cooperative hunting among predators. In both equations in (1), this is expressed by the term $(1+\alpha v)uv$, where $\alpha$ is the predator’s attack rate, which also reflects the degree of cooperative behavior, acting as an Allee threshold. Additionally, in the second equation in (1), the parameter $\sigma$ is included as a food conversion coefficient.

–

Predator mortality is modeled as density-dependent, influenced also by internal processes within the predator population such as competition (or potential cannibalism) (with a mortality rate ${\mu }_2$) and the effects of higher-level predators not explicitly included in the model (with a mortality rate ${\mu }_3$). This is supported by studies in [55,57,58], and is included in the second equation of (1) through the term ${\mu }_1+{\mu }_{2}v+{\mu }_3{v}^2$, as ${\mu }_1$ is the natural death rate of the predator. It is worth noting that this assumption can further limit the predator population size.

–

Like [55], the diffusion coefficients for prey and predator are assumed to differ, as the parameter k in the second equation of (1) accounts for the ratio between the two diffusion constants.

–

The time-fractional derivatives in Equation (1) capture long-term dependencies and correlations between predator and prey populations. These are particularly relevant when the interactions between the species are influenced by events from the more distant past. Moreover, the time-fractional derivatives can also account for memory effects that cannot be established by the standard version of the predator–prey system.

In this study, similar to the research carried out in [54], we shall search for analytical solutions of System (1) assuming that the population waves of the predator and the prey propagate at different speeds. For this purpose, we shall apply the extended SEsM to Equation (1) in two different ways: (1) by introducing fractional traveling wave transformation and using the specter of available integer-order ODEs with known solutions (in Section 3); and (2) by introducing standard traveling wave transformation and using a version of the fractional ODE of Riccati with known solutions (in Section 4). Regarding Case 1, presented above, some key ideas related to the current study were outlined in [54], which I further develop and expand upon here.

2. Description of Extended SEsM Algorithm

In this section, I shall present the extended SEsM algorithm, including several variants for obtaining exact solutions of systems of FNPDEs, like Equation (1). As will be shown below, in contrast to the original SEsM algorithm presented in [51,52,53], I exchange its first and second steps, as extensions in the methodology are made to both steps of the original SEsM algorithm.

Let us consider a system of two FNPDEs with both time- and space-fractional derivatives:

$$\begin{array}{ccc}\hfill & & D_t^{\varphi }u(x,t)+D_x^{\psi }u(x,t)={\Psi }_1(u(x,t),v(x,t),u_x,v_x,u_{xx},\dots \dots )\hfill \\ & & D_t^{\varphi }v(x,t)+D_x^{\psi }v(x,t)={\Psi }_2(u(x,t),v(x,t),u_x,v_x,v_{xx},\dots \dots ),0<\varphi ,\psi \le 1\hfill \end{array}$$

(2)

where $u(x,t)$ and $v(x,t)$ are unknown functions to be solved, x and t are independent variables with respect to space and time, ${\Psi }_1$ and ${\Psi }_2$ are nonlinear functions of $u(x,t)$ and $v(x,t)$ and their spatial derivatives, and $D_t^{\varphi }$ and $D_x^{\psi }$ are the time-fractional derivative of order $\varphi$ and the space-fractional derivative of order $\psi$, respectively.

The extended SEsM algorithm, which is adopted for extracting exact analytical solutions of systems of FNPDEs like (2), involves the following basic steps:

• Construction of the solutions of Equation (2). Apart from the traditional construction of solutions of Equation (2), where the solution of one simple (auxiliary) equation with one independent variable is used (e.g., [41,42,43,44,45]), the SEsM also offers the following possible variant for constructing these solutions:

  $$u\left({\xi }_1\right)=F_1\left({\xi }_1\right),v\left({\xi }_2\right)=F_2\left({\xi }_2\right),$$

  (3)

  where

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

  (4)

  as $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ are the solutions of two simple equations whose solutions are known.
• Selection of the traveling wave type transformation.
  • Variant 1: Use a fractional transformation. The choice of the explicit form of the fractional traveling wave transformation depends on how the fractional derivatives in Equation (2) are defined at the beginning. Below, the most used fractional traveling wave transformations are selected:

    –

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

    –

    Fractional complex transform: ${\xi }_1={\kappa }_1{\displaystyle \frac{x^{\beta }}{\Gamma (\beta +1)}}+{\omega }_1{\displaystyle \frac{t^{\alpha }}{\Gamma (\alpha +1)}},{\xi }_2={\kappa }_2{\displaystyle \frac{x^{\beta }}{\Gamma (\beta +1)}}+{\omega }_2{\displaystyle \frac{t^{\alpha }}{\Gamma (\alpha +1)}}\dots$, defined for modified Riemann–Liouville fractional derivatives [47], which can applied for for Caputo fractional derivatives and other fractional derivatives in studying FNPDEs [59].

    In the both cases, the studied FNPDEs are reduced to integer-order nonlinear ODEs.
  • Variant 2: Use a standard traveling wave transformation. In this case, the introduction of traveling wave ansatz ${\xi }_1={\kappa }_{1}x+{\omega }_{1}t,{\xi }_2={\kappa }_{2}x+{\omega }_{2}t,\dots$ in Equation (2) reduces the studied FNPDEs to fractional nonlinear ODEs.
• Selection of the forms of the simple equations.
  • For Variant 1 of Step 2: The generalized form of the integer-order simple equations used is as follows:

    $${\left({\displaystyle \frac{d^{k_1}{f}_1}{d{\xi }_1^{k_1}}}\right)}^{l_1}=\sum _{j_1=0}^{m_1}{c}_{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}{d}_{j_2}{f}_2^{j_2},$$

    (5)

    where $k_{1,2}$ are the orders of derivatives of $f_1$ and $f_2$, $l_{1,2}$ are the degrees of derivatives in the defining ODEs, and $m_{1,2}$ are the highest degrees of the polynomials of $f_1$ and $f_2$ in the defining ODEs.
    By fixing $l_1$ and $l_2$ and $m_1$ and $m_2$ (at $k_1=k_2=1$) in Equation (5), different types of integer-order ODEs can be used as simple equations, such as the following:

    –

    ODEs of the first order with known analytical solutions (for example, an ODE of Riccati, an ODE of Bernoulli, an ODE of Abel of the first kind, an ODE of the tanh-function, etc.);

    –

    ODEs of the second order with known analytical solutions (for example, elliptic equations of Jaccobi and Weiershtrass and their sub-variants, an ODE of Abel of the 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}{c}_{j_1}{f}_1^{j_1},D_{{\xi }_2}^{\varphi }=\sum _{j_2=0}^{m_2}{d}_{j_2}{f}_2^{j_2}.$$

    (6)
    By fixing $m_1$ and $m_2$ in Equation (6), different types of fractional ODEs can be used as simple equations, such as a fractional ODE of Riccati, a fractional ODE of Bernoulli, or other polynomial or elliptic form equations, whose solutions are available in the literature.
    Note: The correct choice of simple equations used depends on the type of waves that are expected to be realistic for the real-world process being modeled by the corresponding FNPDEs. This is because different types of ODEs include in their solutions different types of special functions, such as hyperbolic, elliptic, or trigonometric functions. Each of these functions describes different wave dynamics of the studied FNPDEs, such as the propagation of solitary waves, kink-type waves, periodic waves, or wave packets. Thus, selecting the right simple equations is essential for accurately describing the system’s behavior regarding the studied real-world processes.
• Derivation of the balance equations and the system of algebraic equations. The fixation of the explicit form of constructed solutions of Equation (2) presented in Step 1 of the SEsM algorithm depends on the balance equations derived. Substitutions of the explicit analytical form of the selected variant equations presented in Step 2 and Step 3 of the SEsM algorithm leads to 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 (2) 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 the formation of a system of nonlinear algebraic equations for each variant chosen according Steps 2 and 3 of the SEsM algorithm.
• Derivation of the analytical solutions. Any non-trivial solution of the algebraic system mentioned above leads to a solution of the studied FNPDEs by replacing the specific expressions of the coefficients in Equations (3) and (4), as well as by changing the traveling wave coordinates chosen from the variants given in Step 2. For simplicity, these solutions are expressed through special functions. For Variant 1 of Step 3, these special functions are $V_{c_0,c_1,\dots ,c_{j_1}}({\xi }_1;k_1,l_1,m_1)$ and $V_{d_0,d_1,\dots ,d_{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 (5)). For Variant 2 of Step 3, the special functions are $V_{c_0,c_1,\dots ,c_{j_1}}({\xi }_1,\varphi ;m_1)$ and $V_{d_0,d_1,\dots ,d_{j_2}}({\xi }_2,\varphi ;m_2)$, whose exact forms are determined by the type of fractional simple equations used (for reference, see Equation (6)).

Following the SEsM algorithm, in the next two sections, I will derive numerous exact solutions of (1). Similar to [60], the obtained analytical results will be stated in terms of lemmas and propositions.

3. Exact Solutions of Time-Fractional Diffusive Predator–Prey System Incorporating Allee Effect Using Fractional Transformation

I present the general solution of Equation (1) as follows:

$$\begin{array}{c}\hfill u\left({\xi }_1\right)=\sum _{i_1=1}^{n_1}{a}_{i_1}{\left[f_1\left({\xi }_1\right)\right]}^{i_1},\\ \hfill v\left({\xi }_2\right)=\sum _{i_2=1}^{n_2}{b}_{i_2}{\left[f_2\left({\xi }_2\right)\right]}^{i_2},\end{array}$$

(7)

In view of the biological nature of the studied time-fractional predator–prey model, we choose to present the memory effects occurring in System (1) in the sense of a Caputo derivative [59]. Thus, the fractional transformations take the following form:

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

(8)

We choose the simple equations used to be nonlinear ODEs of the first order:

$${\displaystyle \frac{d{f}_1}{d{\xi }_1}}=\sum _{j_1=0}^{m_1}{c}_{j_1}{\left[f_1\left({\xi }_1\right)\right]}^{j_1},{\displaystyle \frac{d{f}_2}{d{\xi }_2}}=\sum _{j_2=0}^{m_2}{d}_{j_2}{\left[f_2\left({\xi }_2\right)\right]}^{j_2}$$

(9)

For simplicity, below, we shall present these solutions by the special functions $V_{c_0,c_1,\dots c_{j_1}}$$({\xi }_1;1,1,m_1)$ and $V_{d_0,d_1,\dots d_{j_2}}({\xi }_2;1,1,m_2)$, which are written in the context of Equation (5).

Due to the nature of the studied model, we assume also that the two simple equations used have an identical form. The balance equations are $n_1=m_1-1,n_2=m_2-1$. Below, we shall present various solutions of Equation (1) depending on the numerical values of $m_1$ and $m_2$ in Equation (9).

3.1. Case 1: When $m_1=3$ and $m_2=3$ in Equation (9)

Let us consider the case where $m_1=3$ and $m_2=3$ in Equation (9). Then, the simple equations used are reduced to

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

(10)

3.1.1. Variant 1: When $c_0\ne 0,c_1\ne 0,c_2\ne 0,c_3\ne 0$ and $d_0\ne 0,d_1\ne 0,d_2\ne 0,d_3\ne 0$ in Equation (10)

Let us consider the case where $c_0\ne 0,c_1\ne 0,c_2\ne 0,c_3\ne 0$ and $d_0\ne 0,d_1\ne 0,d_2\ne 0$, $d_3\ne 0$ in Equation (10). For this case, we present the general solution of Equation (1) by special functions in the following manner:

$$\begin{array}{c}\hfill u\left({\xi }_1\right)=a_0+a_1{V}_{c_0,c_1,c_2,c_3}({\xi }_1;1,1,3)+a_2{V}_{c_0,c_1,c_2,c_3}^2({\xi }_1;1,1,3),\\ \hfill v\left({\xi }_2\right)=b_0+b_1{V}_{d_0,d_1,d_2,d_3}({\xi }_2;1,1,3)+b_2{V}_{d_0,d_1,d_2,d_3}^2({\xi }_2;1,1,3),\end{array}$$

(11)

where the special functions

$$\begin{array}{c}\hfill V_{c_0,c_1,c_2,c_3}({\xi }_1;1,1,3)={\displaystyle \frac{exp\left[\left(c_1-\frac{c_2^2}{3c_3}\right)({\xi }_1+{\xi }_{01})\right]}{\sqrt{C_1^{*}-c_{3}exp\left[2\left(c_1-\frac{c_2^2}{3c_3}\right)({\xi }_1+{\xi }_{01})\right]}}}-{\displaystyle \frac{c_2}{3c_3}},\end{array}$$

(12)

$$\begin{array}{c}\hfill V_{d_0,d_1,d_2,d_3}({\xi }_2;1,1,3)={\displaystyle \frac{exp\left[\left(d_1-\frac{d_2^2}{3d_3}\right)({\xi }_2+{\xi }_{02})\right]}{\sqrt{C_2^{*}-d_{3}exp\left[2\left(d_1-\frac{d_2^2}{3d_3}\right)({\xi }_2+{\xi }_{02})\right]}}}-{\displaystyle \frac{d_2}{3d_3}},\end{array}$$

are solutions of ODEs of Abel of the first kind for the particular case $c_0={\displaystyle \frac{c_2}{3c_3}}(c_1-{\displaystyle \frac{2c_2^2}{9c_3}}),d_0={\displaystyle \frac{d_2}{3d_3}}(d_1-{\displaystyle \frac{2d_2^2}{9d_3}})$ [39,61]. In Equation (12), $C_1^{*}$ and $C_2^{*}$ are constants of integration. The Abel ODEs are as follows:

$${\displaystyle \frac{d{f}_1}{d{\xi }_1}}=c_0+c_1{f}_1+c_2{f}_1^2+c_3{f}_1^3,{\displaystyle \frac{d{f}_2}{d{\xi }_2}}=d_0+d_1{f}_2+d_2{f}_2^2+d_3{f}_2^3$$

(13)

Lemma 1. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (13). The application of SEsM with simple equations of kind (13) reduces (1) to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}\hfill & & 8a_2{c_3}^2-{a_2}^3=0\hfill \\ & & 3a_1{c_3}^2-3a_1{a_2}^2+14a_2{c}_3{c}_2=0\hfill \\ & & 12a_2{c}_3{c}_1-3{a_1}^2{a}_2+\beta {a_2}^2+5c_3{a}_1{c}_2+6a_2{c_2}^2+{a_2}^2-3a_0{a_2}^2+2{\omega }_1{a}_2{c}_3=0\hfill \\ & & 2\beta a_1{a}_2-{a_1}^3+10c_2{a}_2{c}_1+2a_1{a}_2+2{\omega }_1{a}_2{c}_2+10a_2{c}_3{c}_0+{\omega }_1{a}_1{c}_3\hfill \\ & & -6a_0{a}_1{a}_2+2a_1{c_2}^2+4c_3{a}_1{c}_1=0\hfill \\ & & a_2\alpha {{\mathit{b}}_{\mathit{2}}}^2-\sigma a_2\alpha {{\mathit{b}}_{\mathit{1}}}^2=0\hfill \\ & & -2a_2\alpha b_0{b}_2-2\sigma a_2\alpha b_0{b}_2-\sigma a_2{b}_2-a_2\alpha {b_1}^2-a_2{b}_2-\sigma a_2\alpha {b_1}^2=0\hfill \\ & & -\sigma a_2{b}_1-a_2{b}_1-2a_2\alpha b_0{b}_1-2\sigma a_2\alpha b_0{b}_1=0\hfill \\ & & 2\beta a_0{a}_2-\sigma a_2\alpha {b_0}^2-a_2\alpha {b_0}^2-\sigma a_2{b}_0+4a_2{c_1}^2-\beta a_2+2a_0{a}_2+\beta {a_1}^2+3c_2{a}_1{c}_1\hfill \\ & & +8c_2{a}_2{c}_0+3c_0{a}_1{c}_3+{\omega }_1{a}_1{c}_2+2{\omega }_1{a}_2{c}_1-3{a_0}^2{a}_2-3a_0{a_1}^2-a_2{b}_0+{a_1}^2=0\hfill \\ & & -a_1\alpha {b_1}^2-2\sigma a_1\alpha b_0{b}_2-a_1{b}_2-2a_1\alpha b_0{b}_2-\sigma a_1\alpha {b_1}^2-\sigma a_1{b}_2=0\hfill \\ & & -\beta a_1-a_1\alpha {b_0}^2-\sigma a_1{b}_0+a_1{c_1}^2+2c_0{a}_1{c}_2+2a_0{a}_1-3{a_0}^2{a}_1-a_1{b}_0\hfill \\ & & +2\beta a_0{a}_1+6c_1{a}_2{c}_0+{\omega }_1{a}_1{c}_1+2{\omega }_1{a}_2{c}_0-\sigma a_1\alpha {b_0}^2=0\hfill \\ & & {\mu }_3{b_2}^3-8\kappa b_2{d_3}^2=0\hfill \\ & & -14\kappa b_2{d}_3{d}_2+3{\mu }_3{b}_1{b_2}^2-3\kappa b_1{d_3}^2=0\hfill \\ & & -a_0\alpha {b_2}^2+{\mu }_2{b_2}^2-6\kappa b_2{d_2}^2-5\kappa d_3{b}_1{d}_2-2{\omega }_2{b}_2{d}_3-12\kappa b_2{d}_3{d}_1\hfill \\ & & +3{\mu }_3{b_1}^2{b}_2+3{\mu }_3{b}_0{b_2}^2-\sigma a_0\alpha {b_2}^2=0\hfill \\ & & -2{\omega }_2{b}_2{d}_2-4\kappa d_3{b}_1{d}_1-2\kappa b_1{d_2}^2-{\omega }_2{b}_1{d}_3-2\sigma a_0\alpha b_1{b}_2+2{\mu }_2{b}_1{b}_2\hfill \\ & & -2a_0\alpha b_1{b}_2-10\kappa d_2{b}_2{d}_1+{\mu }_3{b_1}^3+6{\mu }_3{b}_0{b}_1{b}_2-10\kappa b_2{d}_3{d}_0=0\hfill \\ & & -3\kappa d_0{b}_1{d}_3+{\mu }_1{b}_2-3\kappa d_2{b}_1{d}_1-2{\omega }_2{b}_2{d}_1+3{\mu }_3{b}_0{b_1}^2\hfill \\ & & -2\sigma a_0\alpha b_0{b}_2+{\mu }_2{b_1}^2-8\kappa d_2{b}_2{d}_0-2a_0\alpha b_0{b}_2\hfill \\ & & -a_0\alpha {b_1}^2-\sigma a_0\alpha {b_1}^2-\sigma a_0{b}_2-a_0{b}_2+3{\mu }_3{b_0}^2{b}_2-4\kappa b_2{d_1}^2+2{\mu }_2{b}_0{b}_2-{\omega }_2{b}_1{d}_2=0\hfill \\ & & -\sigma a_0{b}_1-2a_0\alpha b_0{b}_1+3{\mu }_3{b_0}^2{b}_1-2\sigma a_0\alpha b_0{b}_1-2{\omega }_2{b}_2{d}_0-a_0{b}_1\hfill \\ & & -2\kappa d_0{b}_1{d}_2-\kappa b_1{d_1}^2+{\mu }_1{b}_1+2{\mu }_2{b}_0{b}_1-6\kappa d_1{b}_2{d}_0-{\omega }_2{b}_1{d}_1=0\hfill \\ & & {\mu }_1{b}_0+{\mu }_2{b_0}^2+{\mu }_3{b_0}^3-2\kappa b_2{d_0}^2+2a_2{c_0}^2-\beta a_0+c_0{a}_1{c}_1+{\omega }_1{a}_1{c}_0-a_0{b}_0\hfill \\ & & -\sigma a_0\alpha {b_0}^2+{a_0}^2+\beta {a_0}^2-{a_0}^3-\kappa d_0{b}_1{d}_1-a_0\alpha {b_0}^2-{\omega }_2{b}_1{d}_0-\sigma a_0{b}_0=0\hfill \end{array}$$

(14)

Proof. 

The application of the extended SEsM to (1) means that we have to substitute (7) for the case $n_1=n_2=2$ together with (8) and Equation (13) into Equation (1). Then, we have to equate to 0 the coefficients of the obtained polynomials of $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ in (1). This leads us to the algebraic system presented in (14). □

Proposition 1. 

An exact solution of Equation (1) is

$$\begin{array}{ccc}\hfill & & u^{\left(1\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{18c_3}}\left({\Omega }_1-2\sqrt{2}{c_2}^2+3{\omega }_1\sqrt{2}{c}_3+18\sqrt{2}{c}_3{c}_1+9c_3\right)\hfill \\ & & +{\displaystyle \frac{4}{3}}\sqrt{2}{c}_2{V}_{c_0,c_1,c_2,c_3}({\xi }_1;1,1,3)+2\sqrt{2}{c}_3{V}_{c_0,c_1,c_2,c_3}^2({\xi }_1;1,1,3)\hfill \\ & & v^{\left(1\right)}\left({\xi }_2\right)={\displaystyle \frac{1}{18}}{\displaystyle \frac{-9{\mu }_3{d}_3+8\alpha \sqrt{2}{d_2}^2\sqrt{\kappa {\mu }_3}}{\alpha {\mu }_3{d}_3}}\hfill \\ & & -{\displaystyle \frac{4}{3}}{\displaystyle \frac{\kappa d_2\sqrt{2}}{\sqrt{\kappa {\mu }_3}}}{V}_{d_0,d_1,d_2,d_3}({\xi }_2;1,1,3)-2{\displaystyle \frac{\kappa d_3\sqrt{2}}{\sqrt{\kappa {\mu }_3}}}{V}_{d_0,d_1,d_2,d_3}^2({\xi }_2;1,1,3)\hfill \end{array}$$

(15)

where

$$\begin{array}{c}\hfill V_{c_0,c_1,c_2,c_3}({\xi }_1;1,1,3)={\displaystyle \frac{exp\left[\left(c_1-\frac{c_2^2}{3c_3}\right)({\xi }_1+{\xi }_{01})\right]}{\sqrt{C_1^{*}-c_{3}exp\left[2\left(c_1-\frac{c_2^2}{3c_3}\right)({\xi }_1+{\xi }_{01})\right]}}}-{\displaystyle \frac{c_2}{3c_3}},\\ \hfill V_{d_0,d_1,d_2,d_3}({\xi }_2;1,1,3)={\displaystyle \frac{exp\left[\left(d_1-\frac{d_2^2}{3d_3}\right)({\xi }_2+{\xi }_{02})\right]}{\sqrt{C_2^{*}-d_{3}exp\left[2\left(d_1-\frac{d_2^2}{3d_3}\right)({\xi }_2+{\xi }_{02})\right]}}}-{\displaystyle \frac{d_2}{3d_3}},\end{array}$$

for the particular case $c_0={\displaystyle \frac{c_2}{3c_3}}(c_1-{\displaystyle \frac{2c_2^2}{9c_3}}),d_0={\displaystyle \frac{d_2}{3d_3}}(d_1-{\displaystyle \frac{2d_2^2}{9d_3}})$ and ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}},{\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (13). We apply the extended SEsM to (1) on the basis of substitutions of (7) in the case $n_1=n_2=2$, (8) and (13) into (1). As a result, we obtain the system of nonlinear algebraic equations presented in (14). A non–trivial solution of this system is as follows:

$$\begin{array}{ccc}\hfill & & a_0={\displaystyle \frac{1}{18c_3}}\left({\Omega }_1-2\sqrt{2}{c_2}^2+3{\omega }_1\sqrt{2}{c}_3+18\sqrt{2}{c}_3{c}_1+9c_3\right),a_1={\displaystyle \frac{4}{3}}\sqrt{2}{c}_2,a_2=2\sqrt{2}{c}_3,\hfill \\ & & b_0={\displaystyle \frac{1}{18}}{\displaystyle \frac{-9{\mu }_3{d}_3+8\alpha \sqrt{2}{d_2}^2\sqrt{\kappa {\mu }_3}}{\alpha {\mu }_3{d}_3}},b_1=-{\displaystyle \frac{4}{3}}{\displaystyle \frac{\kappa d_2\sqrt{2}}{\sqrt{\kappa {\mu }_3}}},b_2=-2{\displaystyle \frac{\kappa d_3\sqrt{2}}{\sqrt{\kappa {\mu }_3,}}},\sigma ={\displaystyle \frac{b_2^2}{b_1^2}},\hfill \\ & & {\mu }_1={\displaystyle \frac{1}{36{\alpha }^2\sqrt{\kappa {\mu }_3}{d_3}^2}},(27\sqrt{\kappa {\mu }_3}{d_3}^2{\mu }_3+16{\alpha }^2\sqrt{\kappa {\mu }_3}\kappa {d_2}^4+144{\alpha }^2\sqrt{\kappa {\mu }_3}{d_3}^2\kappa {d_1}^2\hfill \\ & & +72{\alpha }^2\sqrt{\kappa {\mu }_3}{d_3}^2{\omega }_2{d}_1+192{\alpha }^2\sqrt{\kappa {\mu }_3}{d}_3\kappa {d_2}^2{d}_1-18\alpha {\mu }_3{d_3}^2\sqrt{2}{\omega }_2-108\alpha \kappa {\mu }_3{d_3}^2\sqrt{2}{d}_1\hfill \\ & & +24{\alpha }^2\sqrt{\kappa {\mu }_3}{d}_3{d_2}^2{\omega }_2-36\alpha \kappa {\mu }_3{d}_3\sqrt{2}{d_2}^2),\beta ={\displaystyle \frac{1}{2}}{\displaystyle \frac{c_3+\sqrt{{c_3}^2\left(-3+24{c_1}^2+2{{\omega }_1}^2\right)}}{c_3}}\hfill \\ & & {\mu }_2={\displaystyle \frac{1}{4}}{\displaystyle \frac{\left(-12\alpha \kappa d_3\sqrt{\kappa {\mu }_3}{d}_1-4\alpha \kappa \sqrt{\kappa {\mu }_3}{d_2}^2-2\alpha d_3\sqrt{\kappa {\mu }_3}{\omega }_2+3\kappa d_3{\mu }_3\sqrt{2}\right)\sqrt{2}}{\alpha \kappa d_3}}\hfill \end{array}$$

(16)

where

$${\Omega }_1=\sqrt{-27{c_3}^2+18{{\omega }_1}^2{c_3}^2-144{c_2}^2{c}_1{c}_3+216{c_3}^2{c_1}^2+24{c_2}^4}$$

The substitution of $a_0,a_1$, and $a_2$ and $b_0,b_1$, and $b_2$ from (16) into (11), together with (12), leads to the solution (15) of (1). □

3.1.2. Variant 2: When $c_0=0,c_1\ne 0,c_2=0,c_3\ne 0$ and $d_0=0,d_1\ne 0,d_2=0,d_3\ne 0$ in Equation (10)

For this case, we present the general solution of Equation (1) by special functions in the following manner:

$$\begin{array}{c}\hfill u\left({\xi }_1\right)=a_0+a_1{V}_{0,c_1,0,c_3}({\xi }_1;1,1,3)+a_2{V}_{0,c_1,0,c_3}^2({\xi }_1;1,1,3),\\ \hfill v\left({\xi }_2\right)=b_0+b_1{V}_{0,d_1,0,d_3}({\xi }_2;1,1,3)+b_2{V}_{0,d_1,0,d_3}^2({\xi }_2;1,1,3),\end{array}$$

(17)

where the special functions $V_{0,c_1,0,c_3}({\xi }_1;1,1,3)$ and $V_{0,d_1,0,d_3}({\xi }_2;1,1,3)$ are solutions of a particular variants of a Bernoulli ODE [39,50].

$$V_{0,c_1,0,c_3}^{\left(1\right)}({\xi }_1;1,1,3)=\sqrt{{\displaystyle \frac{c_{1}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}{1-c_{3}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}}}$$

(18)

for $c_1>0,c_3<0$,

$$V_{0,d_1,0,d_3}^{\left(1\right)}({\xi }_2;1,1,3)=\sqrt{{\displaystyle \frac{d_{1}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}{1-d_{3}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}}}$$

(19)

for $d_1>0,d_3<0$,

$$V_{0,c_1,0,c_3}^{\left(2\right)}({\xi }_1;1,1,3)=-\sqrt{{\displaystyle \frac{c_{1}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}{1+c_{3}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}}}$$

(20)

for $c_1<0,c_3>0$, and

$$V_{0,d_1,0,d_3}^{\left(2\right)}({\xi }_2;1,1,3)=-\sqrt{{\displaystyle \frac{d_{1}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}{1+d_{3}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}}}$$

(21)

for $d_1<0,d_3>0$. The particular variants of the Bernoulli ODEs are as follows:

$${\displaystyle \frac{d{f}_1}{d{\xi }_1}}=c_1{f}_1+c_3{f}_1^3,{\displaystyle \frac{d{f}_2}{d{\xi }_2}}=d_1{f}_2+d_3{f}_2^3$$

(22)

The combination of different solutions of the ODEs of Bernoulli, depending on the numerical values of their coefficients, leads to obtaining several types of analytical solutions of Equation (1).

Lemma 2. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (22). The application of the extended SEsM with simple equations of kind (22) reduces (1) to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}\hfill & & 8a_2{c_3}^2-{a_2}^3=0\hfill \\ & & -3a_1{a_2}^2+3a_1{c_3}^2=0\hfill \\ & & -{a_1}^3+2a_1{a}_2+{\omega }_1{a}_1{c}_3+4a_1{c}_1{c}_3+2\beta a_1{a}_2-6a_0{a}_1{a}_2=0\hfill \\ & & -a_2\alpha {b_2}^2-\sigma a_2\alpha {b_1}^2=0\hfill \\ & & -2\sigma a_2\alpha b_0{b}_2-\sigma a_2\alpha {b_1}^2-a_2{b}_2-2a_2\alpha b_0{b}_2-\sigma a_2{b}_2-a_2\alpha {b_1}^2=0\hfill \\ & & -a_2{b}_1-2\sigma a_2\alpha b_0{b}_1-\sigma a_2{b}_1-2a_2\alpha b_0{b}_1=0\hfill \\ & & 4a_2{c_1}^2-3{a_0}^2{a}_2-a_2{b}_0-3a_0{a_1}^2+\beta {a_1}^2+2a_0{a}_2-\beta a_2+2{\omega }_1{a}_2{c}_1\hfill \\ & & +2\beta a_0{a}_2-\sigma a_2\alpha {b_0}^2-a_2\alpha {b_0}^2-\sigma a_2{b}_0+{a_1}^2=0\hfill \\ & & -\sigma a_1\alpha {b_0}^2+2\beta a_0{a}_1-a_1\alpha {b_0}^2-a_1{b}_0+{\omega }_1{a}_1{c}_1-\beta a_1+2a_0{a}_1-\sigma a_1{b}_0+a_1{c_1}^2-3{a_0}^2{a}_1=0\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& & {\mu }_3{b_2}^3-8\kappa b_2{d_3}^2=0\hfill \\ & & 3{\mu }_3{b}_1{b_2}^2-3\kappa b_1{d_3}^2=0\hfill \\ & & -\sigma a_0\alpha {b_2}^2+{\mu }_2{b_2}^2+3{\mu }_3{b_1}^2{b}_2+3{\mu }_3{b}_0{b_2}^2-a_0\alpha {b_2}^2-12\kappa b_2{d}_1{d}_3-2{\omega }_2{b}_2{d}_3=0\hfill \\ & & 2{\mu }_2{b}_1{b}_2+{\mu }_3{b_1}^3-{\omega }_2{b}_1{d}_3-4\kappa b_1{d}_1{d}_3-2a_0\alpha b_1{b}_2+6{\mu }_3{b}_0{b}_1{b}_2-2\sigma a_0\alpha b_1{b}_2=0\hfill \\ & & -2{\omega }_2{b}_2{d}_1-\sigma a_0\alpha {b_1}^2-2a_0\alpha b_0{b}_2+3{\mu }_3{b_0}^2{b}_2+{\mu }_2{b_1}^2-2\sigma a_0\alpha b_0{b}_2\hfill \\ & & -4\kappa b_2{d_1}^2-a_0{b}_2+{\mu }_1{b}_2+2{\mu }_2{b}_0{b}_2+3{\mu }_3{b}_0{b_1}^2-\sigma a_0{b}_2-a_0\alpha {b_1}^2=0\hfill \\ & & -a_0{b}_1+3{\mu }_3{b_0}^2{b}_1+2{\mu }_2{b}_0{b}_1+{\mu }_1{b}_1-\sigma a_0{b}_1-2a_0\alpha b_0{b}_1\hfill \\ & & -\kappa b_1{d_1}^2-{\omega }_2{b}_1{d}_1-2\sigma a_0\alpha b_0{b}_1=0\hfill \\ & & \beta a_0+{a_0}^2+\beta {a_0}^2-a_0{b}_0+{\mu }_1{b}_0+{\mu }_2{b_0}^2+{\mu }_3{b_0}^3-\sigma a_0\alpha {b_0}^2-{a_0}^3-a_0\alpha {b_0}^2-\sigma a_0{b}_0=0\hfill \end{array}$$

Proof. 

The application of the extended SEsM to (1) means that we have to substitute (7) for the case $n_1=n_2=2$, together with (8) and (22) into (1). Then, we have to equate to 0 the coefficients of the obtained polynomials of $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ in (1). This leads us to the algebraic system presented in (23). □

Proposition 2. 

The exact solutions of Equation (1) are as follows:

$$\begin{array}{ccc}\hfill & & u^{\left(2\right)}\left({\xi }_1\right)=-{\displaystyle \frac{1}{15}}{\omega }_1+{\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{90}}{\Omega }_1+{\displaystyle \frac{1}{15}}\sqrt{15}\sqrt{{\omega }_1{c}_3}{V}_{0,c_1,0,c_3}^{\left(1\right)}({\xi }_1;1,1,3)-c_3{V^{\left(1\right)}}_{0,c_1,0,c_3}^2({\xi }_1;1,1,3),\hfill \\ & & v^{\left(2\right)}\left({\xi }_2\right)={\displaystyle \frac{1}{450}}{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }\sqrt{5}\sqrt{-{\mu }_3{\Omega }_2{\Omega }_3}+150{{\mu }_3}^2\kappa {{\omega }_2}^2}{{{\mu }_3}^2\kappa {\Omega }_2}}\hfill \\ & & +2{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_3}{{\mu }_3}}{V^{\left(1\right)}}_{0,d_1,0,d_3}^2({\xi }_2;1,1,3),\hfill \end{array}$$

(24)

where

$$V_{0,c_1,0,c_3}^{\left(1\right)}({\xi }_1;1,1,3)=\sqrt{{\displaystyle \frac{c_{1}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}{1-c_{3}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}}}$$

for $c_1>0,c_3<0$,

$$V_{0,d_1,0,d_3}^{\left(1\right)}({\xi }_2;1,1,3)=\sqrt{{\displaystyle \frac{d_{1}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}{1-d_{3}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}}}$$

for $d_1>0,d_3<0$, with ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}}$ and ${\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

$$\begin{array}{ccc}\hfill & & u^{\left(3\right)}\left({\xi }_1\right)=-{\displaystyle \frac{1}{15}}{\omega }_1+{\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{90}}{\Omega }_1+{\displaystyle \frac{1}{15}}\sqrt{15}\sqrt{{\omega }_1{c}_3}{V}_{0,c_1,0,c_3}^{\left(2\right)}({\xi }_1;1,1,3)-c_3{V^{\left(2\right)}}_{0,c_1,0,c_3}^2({\xi }_1;1,1,3),\hfill \\ & & v^{\left(3\right)}\left({\xi }_2\right)={\displaystyle \frac{1}{450}}{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }\sqrt{5}\sqrt{-{\mu }_3{\Omega }_2{\Omega }_3}+150{{\mu }_3}^2\kappa {{\omega }_2}^2}{{{\mu }_3}^2\kappa {\Omega }_2}}\hfill \\ & & +2{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_3}{{\mu }_3}}{V^{\left(2\right)}}_{0,d_1,0,d_3}^2({\xi }_2;1,1,3),\hfill \end{array}$$

(25)

where

$$V_{0,c_1,0,c_3}^{\left(2\right)}({\xi }_1;1,1,3)=-\sqrt{{\displaystyle \frac{c_{1}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}{1+c_{3}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}}}$$

for $c_1<0,c_3>0$,

$$V_{0,d_1,0,d_3}^{\left(2\right)}({\xi }_2;1,1,3)=-\sqrt{{\displaystyle \frac{d_{1}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}{1+d_{3}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}}}$$

for $d_1<0,d_3>0$, where ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}}$ and ${\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

$$\begin{array}{ccc}\hfill & & u^{\left(4\right)}\left({\xi }_1\right)=-{\displaystyle \frac{1}{15}}{\omega }_1+{\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{90}}{\Omega }_1+{\displaystyle \frac{1}{15}}\sqrt{15}\sqrt{{\omega }_1{c}_3}{V}_{0,c_1,0,c_3}^{\left(1\right)}({\xi }_1;1,1,3)-c_3{V^{\left(1\right)}}_{0,c_1,0,c_3}^2({\xi }_1;1,1,3),\hfill \\ & & v^{\left(4\right)}\left({\xi }_2\right)={\displaystyle \frac{1}{450}}{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }\sqrt{5}\sqrt{-{\mu }_3{\Omega }_2{\Omega }_3}+150{{\mu }_3}^2\kappa {{\omega }_2}^2}{{{\mu }_3}^2\kappa {\Omega }_2}}\hfill \\ & & +2{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_3}{{\mu }_3}}{V^{\left(2\right)}}_{0,d_1,0,d_3}^2({\xi }_2;1,1,3),\hfill \end{array}$$

(26)

where

$$V_{0,c_1,0,c_3}^{\left(1\right)}({\xi }_1;1,1,3)=\sqrt{{\displaystyle \frac{c_{1}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}{1-c_{3}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}}}$$

for $c_1>0,c_3<0$,

$$V_{0,d_1,0,d_3}^{\left(2\right)}({\xi }_2;1,1,3)=-\sqrt{{\displaystyle \frac{d_{1}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}{1+d_{3}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}}}$$

for $d_1<0,d_3>0$, ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}}$, and ${\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

$$\begin{array}{ccc}\hfill & & u^{\left(5\right)}\left({\xi }_1\right)=-{\displaystyle \frac{1}{15}}{\omega }_1+{\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{90}}{\Omega }_1+{\displaystyle \frac{1}{15}}\sqrt{15}\sqrt{{\omega }_1{c}_3}{V}_{0,c_1,0,c_3}^{\left(2\right)}({\xi }_1;1,1,3)-c_3{V^{\left(2\right)}}_{0,c_1,0,c_3}^2({\xi }_1;1,1,3),\hfill \\ & & v^{\left(5\right)}\left({\xi }_2\right)={\displaystyle \frac{1}{450}}{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }\sqrt{5}\sqrt{-{\mu }_3{\Omega }_2{\Omega }_3}+150{{\mu }_3}^2\kappa {{\omega }_2}^2}{{{\mu }_3}^2\kappa {\Omega }_2}}\hfill \\ & & +2{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_3}{{\mu }_3}}{V^{\left(1\right)}}_{0,d_1,0,d_3}^2({\xi }_2;1,1,3),\hfill \end{array}$$

(27)

where

$$V_{0,c_1,0,c_3}^{\left(2\right)}({\xi }_1;1,1,3)=-\sqrt{{\displaystyle \frac{c_{1}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}{1+c_{3}exp\left[2c_1({\xi }_1+{\xi }_{01})\right]}}}$$

for $c_1<0,c_3>0$,

$$V_{0,d_1,0,d_3}^{\left(1\right)}({\xi }_2;1,1,3)=\sqrt{{\displaystyle \frac{d_{1}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}{1-d_{3}exp\left[2d_1({\xi }_2+{\xi }_{02})\right]}}}$$

for $d_1>0,d_3<0$, ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}}$ and ${\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (22). We apply the extended SEsM to (1) on the basis of (7) in the case $n_1=n_2=2$, (8) and (22). As a result, we obtain the system of nonlinear algebraic equations presented in (23). A non–trivial solution of this system is

$$\begin{array}{ccc}\hfill & & a_0=-{\displaystyle \frac{1}{15}}{\omega }_1+{\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{90}}{\Omega }_1,a_1={\displaystyle \frac{1}{15}}\sqrt{15}\sqrt{{\omega }_1{c}_3},a_2=-c_3,\hfill \\ & & b_0={\displaystyle \frac{1}{450}}{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }\sqrt{5}\sqrt{-{\mu }_3{\Omega }_2{\Omega }_3}+150{{\mu }_3}^2\kappa {{\omega }_2}^2}{{{\mu }_3}^2\kappa {\Omega }_2}},b_1=0,b_2=2{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_3}{{\mu }_3}}\hfill \end{array}$$

(28)

$$\begin{array}{ccc}& & c_1={\displaystyle \frac{1}{5}}\left(-{\displaystyle \frac{1}{6}}+{\displaystyle \frac{1}{6}}\sqrt{37}\right){\omega }_1,d_1={\displaystyle \frac{1}{450}}{\displaystyle \frac{\sqrt{-5{\mu }_3{\Omega }_2{\Omega }_3}}{{\mu }_3{\Omega }_2\kappa }},\beta ={\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{30}}{\Omega }_1,\hfill \\ & & {\mu }_1={\displaystyle \frac{1}{10125}}{\displaystyle \frac{-2250{{\omega }_2}^3\sqrt{2}{\mu }_3\sqrt{{\mu }_3\kappa }+225{\Omega }_1{{\mu }_3}^2{\kappa }^2+3375{{\mu }_2}^2\kappa {\omega }_2\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3{\Omega }_2\kappa }}\hfill \end{array}$$

where

$$\begin{array}{ccc}& & {\Omega }_1=\sqrt{-675{\alpha }_2+12{{\omega }_1}^2-84{{\omega }_1}^2\alpha \sqrt{37}},{\Omega }_2=\kappa {\mu }_2+{\omega }_2\sqrt{2}\sqrt{{\mu }_3\kappa },\hfill \\ & & {\Omega }_3=450{\Omega }_1{{\mu }_3}^2{\kappa }^2-2250{{\mu }_2}^3{\kappa }^2+14{\Omega }_1{{\omega }_1}^2\sqrt{37}{{\mu }_3}^2{\kappa }^2+3375\kappa {\mu }_2{\mu }_3{{\omega }_2}^2\hfill \end{array}$$

The substitution of $a_0,a_1$, and $a_2$ and $b_0,b_1$, and $b_2$ from (28) into (11), together with (18)–(21), leads to the solutions (24)–(27) of (1). □

3.2. Case 2: When $m_1=2$ and $m_2=2$ in Equation (9)

Let us consider the case where $m_1=2$ and $m_2=2$ in Equation (9). Then, the simple equations used are reduced to

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

(29)

3.2.1. Variant 1: When $c_0\ne 0,c_1\ne 0,c_2\ne 0$ and $d_0\ne 0,d_1\ne 0,d_2\ne 0$ in Equation (29)

For this case, the general solution of Equation (1) can be presented as

$$\begin{array}{c}\hfill u\left({\xi }_1\right)=a_0+a_1{V}_{c_0,c_1,c_2}({\xi }_1;1,1,2),\\ \hfill v\left({\xi }_2\right)=b_0+b_1{V}_{d_0,d_1,d_2}({\xi }_2;1,1,2),\end{array}$$

(30)

where the special functions $V_{c_0,c_1,c_2}({\xi }_1;1,1,2)$ and $V_{d_0,d_1,d_2}({\xi }_2;1,1,2)$ are the general solutions of an ODE of Riccati [39,50,60].

$$\begin{array}{c}\hfill V_{c_0,c_1,c_2}({\xi }_1;1,1,2)={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\mathit{c}}_{\mathit{1}}}{{\mathit{c}}_{\mathit{2}}}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\Delta }_1}{{\mathit{c}}_{\mathit{2}}}}tanh\left[{\displaystyle \frac{{\Delta }_1}{2}}({\xi }_1+{\xi }_{01})\right]\\ \hfill +{\displaystyle \frac{{\mathit{F}}_{\mathit{1}}}{{cosh}^2\left[\frac{{\Delta }_1}{2}({\xi }_1+{\xi }_{01})\right]\left({\mathit{E}}_{\mathit{1}}-2\frac{{\mathit{c}}_{\mathit{2}}{\mathit{F}}_{\mathit{1}}}{{\Delta }_1}tanh\left[\frac{{\Delta }_1}{2}({\xi }_1+{\xi }_{01})\right]\right)}}\\ \hfill V_{d_0,d_1,d_2}({\xi }_2;1,1,2)={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\mathit{d}}_{\mathit{1}}}{{\mathit{d}}_{\mathit{2}}}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\Delta }_2}{{\mathit{d}}_{\mathit{2}}}}tanh\left[{\displaystyle \frac{{\Delta }_2}{2}}({\xi }_2+{\xi }_{02})\right]\\ \hfill +{\displaystyle \frac{{\mathit{F}}_{\mathit{2}}}{{cosh}^2\left[\frac{{\Delta }_2}{2}({\xi }_2+{\xi }_{02})\right]\left({\mathit{E}}_{\mathit{2}}-2\frac{{\mathit{d}}_{\mathit{2}}{\mathit{F}}_{\mathit{2}}}{{\Delta }_2}tanh\left[\frac{{\Delta }_2}{2}({\xi }_2+{\xi }_{02})\right]\right)}}\end{array}$$

(31)

where ${\Delta }_1=c_1^2-4c_0{c}_2$ and ${\Delta }_2=d_1^2-4d_0{d}_2$, and $F_1,F_2,E_1$, and $E_2$ are constants. The Riccati ODEs are

$${\displaystyle \frac{d{f}_1}{d{\xi }_1}}=c_0+c_1{f}_1+c_2{f}_1^2,{\displaystyle \frac{d{f}_2}{d{\xi }_2}}=d_0+d_1{f}_2+d_2{f}_2^2$$

(32)

Lemma 3. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (32). The application of the extended SEsM with simple equations of kind (32) reduces Equation (1) to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}\hfill & & 2a_1{c_2}^2-{a_1}^3=0\hfill \\ & & -3a_0{a_1}^2+{a_1}^2+\beta {a_1}^2+3a_1{c}_1{c}_2+{\omega }_1{a}_1{c}_2=0\hfill \\ & & -\sigma a_1\alpha {b_1}^2-a_1\alpha {b_2}^2=0\hfill \\ & & -a_1{b}_1-\sigma a_1{b}_1-2\sigma a_1\alpha b_0{b}_1-2a_1\alpha b_0{b}_1=0\hfill \\ & & 2\beta a_0{a}_1-\sigma a_1\alpha {b_0}^2-\sigma a_1{b}_0-a_1{b}_0-\beta a_1+2a_0{a}_1+{\omega }_1{a}_1{c}_1\hfill \\ & & -a_1\alpha {b_0}^2+a_1{c_1}^2-3{a_0}^2{a}_1+2a_1{c}_2{c}_0=0\hfill \\ & & {\mu }_3{b_1}^3-2\kappa b_1{d_2}^2=0\hfill \\ & & -a_0\alpha {b_1}^2+3{\mu }_3{b}_0{b_1}^2-\sigma a_0\alpha {b_1}^2-{\omega }_2{b}_1{d}_2-3\kappa b_1{d}_1{d}_2+{\mu }_2{b_1}^2=0\hfill \\ & & -2\kappa b_1{d}_2{d}_0-a_0{b}_1+{\mu }_1{b}_1-\kappa b_1{d_1}^2+3{\mu }_3{b_0}^2{b}_1-{\omega }_2{b}_1{d}_1\hfill \\ & & -\sigma a_0{b}_1-2a_0\alpha b_0{b}_1+2{\mu }_2{b}_0{b}_1-2\sigma a_0\alpha b_0{b}_1=0\hfill \\ & & -a_0\alpha {b_0}^2+a_1{c}_1{c}_0-a_0{b}_0+{\mu }_1{b}_0-\sigma a_0{b}_0+{\mu }_3{b_0}^3-{a_0}^3\hfill \\ & & -\sigma a_0\alpha {b_0}^2+{\omega }_1{a}_1{c}_0+\beta {a_0}^2-{\omega }_2{b}_1{d}_0+{a_0}^2+{\mu }_2{b_0}^2-\kappa b_1{d}_1{d}_0-\beta a_0=0\hfill \end{array}$$

(33)

Proof. 

The application of the extended SEsM to Equation (1) means that we have to substitute (7) for the case $n_1=n_2=1$, together with (8) and (32) into (1). Then, we have to equate to 0 the coefficients of the obtained polynomials of $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ in Equation (1). This leads us to the algebraic system presented in Equation (33). □

Proposition 3. 

An exact solution of Equation (1) is

$$\begin{array}{ccc}\hfill & & u^{\left(6\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}{c}_1+{\displaystyle \frac{1}{6}}{\Omega }_4+\sqrt{2}{c}_2{V}_{c_0,c_1,c_2}({\xi }_1;1,1,2),\hfill \\ & & v^{\left(6\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{6}}{\displaystyle \frac{{\omega }_2\sqrt{2}\sqrt{{\mu }_3\kappa }+2{\mu }_2\kappa +3\kappa \sqrt{2}\sqrt{{\mu }_3\kappa }{d}_1}{{\mu }_3\kappa }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_2}{{\mu }_3}}{V}_{d_0,d_1,d_2}({\xi }_2;1,1,2),\hfill \end{array}$$

(34)

where

$$\begin{array}{c}\hfill V_{c_0,c_1,c_2}({\xi }_1;1,1,2)={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\mathit{c}}_{\mathit{1}}}{{\mathit{c}}_{\mathit{2}}}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\Delta }_1}{{\mathit{c}}_{\mathit{2}}}}tanh\left[{\displaystyle \frac{{\Delta }_1}{2}}({\xi }_1+{\xi }_{01})\right]\\ \hfill +{\displaystyle \frac{{\mathit{F}}_{\mathit{1}}}{{cosh}^2\left[\frac{{\Delta }_1}{2}({\xi }_1+{\xi }_{01})\right]\left({\mathit{E}}_{\mathit{1}}-2\frac{{\mathit{c}}_{\mathit{2}}{\mathit{F}}_{\mathit{1}}}{{\Delta }_1}tanh\left[\frac{{\Delta }_1}{2}({\xi }_1+{\xi }_{01})\right]\right)}}\\ \hfill V_{d_0,d_1,d_2}({\xi }_2;1,1,2)={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\mathit{d}}_{\mathit{1}}}{{\mathit{d}}_{\mathit{2}}}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\Delta }_2}{{\mathit{d}}_{\mathit{2}}}}tanh\left[{\displaystyle \frac{{\Delta }_2}{2}}({\xi }_2+{\xi }_{02})\right]\\ \hfill +{\displaystyle \frac{{\mathit{F}}_{\mathit{2}}}{{cosh}^2\left[\frac{{\Delta }_2}{2}({\xi }_2+{\xi }_{02})\right]\left({\mathit{E}}_{\mathit{2}}-2\frac{{\mathit{d}}_{\mathit{2}}{\mathit{F}}_{\mathit{2}}}{{\Delta }_2}tanh\left[\frac{{\Delta }_2}{2}({\xi }_2+{\xi }_{02})\right]\right)}}\end{array}$$

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

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of Equation (1) are of kind (32). We apply the extended SEsM to Equation (1) on the basis of (7) for the case $n_1=n_2=1$, as well as (8) and (32). As a result, we obtain the system of nonlinear algebraic equations presented in (33). A non–trivial solution of this system is

$$\begin{array}{ccc}\hfill & & a_0={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}{c}_1+{\displaystyle \frac{1}{6}}{\Omega }_4,a_1=\sqrt{2}{c}_2,b_1=-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_2}{{\mu }_3}}\hfill \\ & & b_0=-{\displaystyle \frac{1}{6}}{\displaystyle \frac{{\omega }_2\sqrt{2}\sqrt{{\mu }_3\kappa }+2{\mu }_2\kappa +3\kappa \sqrt{2}\sqrt{{\mu }_3\kappa }{d}_1}{{\mu }_3\kappa }},\beta ={\displaystyle \frac{1+{\Omega }_4}{2}},\sigma ={\displaystyle \frac{b_2^2}{b_1^2}},\hfill \\ & & d_0={\displaystyle \frac{1}{36}}{\displaystyle \frac{-{{\omega }_2}^3\sqrt{2}{\mu }_3\sqrt{{\mu }_3\kappa }+2{{\mu }_2}^3{\kappa }^2-3{\Omega }_4{c_1}^2{{\mu }_3}^2{\kappa }^2-{\Omega }_4{{\omega }_1}^2{{\mu }_3}^2{\kappa }^2+12{\Omega }_4{c}_2{c}_0{{\mu }_3}^2{\kappa }^2}{{\Omega }_5{\kappa }^2{d}_2{\mu }_3}},\hfill \\ & & {\mu }_1={\displaystyle \frac{1}{18}}{\displaystyle \frac{-4{{\mu }_2}^3{\kappa }^2-{{\omega }_1}^3\sqrt{2}{{\mu }_3}^2{\kappa }^2-4{{\omega }_2}^3\sqrt{2}{\mu }_3\sqrt{{\mu }_3\kappa }+3{\Omega }_4-{\Omega }_4{{\omega }_1}^2{{\mu }_3}^2{\kappa }^2}{{\mu }_3\kappa {\Omega }_5}}\hfill \end{array}$$

(35)

where

$${\Omega }_4=\sqrt{-3+2{{\omega }_1}^2\alpha +6{c_1}^2-24c_2{c}_0{\alpha }^2},{\Omega }_5=-{\mu }_2\kappa +{\omega }_2\sqrt{2}\sqrt{{\mu }_3\kappa }$$

Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (35) into (30), together with (21), leads to (34) as the solution to (1). □

3.2.2. Variant 2: When $c_0=0,c_1\ne 0,c_2\ne 0$ and $d_0=0,d_1\ne 0,d_2\ne 0$ in Equation (29)

For this case, the general solution of Equation (1) can be presented as

$$\begin{array}{c}\hfill u\left({\xi }_1\right)=a_0+a_1{V}_{0,c_1,c_2}({\xi }_1;1,1,2),\\ \hfill v\left({\xi }_2\right)=b_0+b_1{V}_{0,d_1,d_2}({\xi }_2;1,1,2),\end{array}$$

(36)

where the special functions $V_{0,c_1,c_2}({\xi }_1;1,1,2)$ and $V_{0,d_1,d_2}({\xi }_2;1,1,2)$ are the solutions of particular variants of a Bernoulli ODE [39,50].

$$V_{0,c_1,c_2}^{\left(1\right)}({\xi }_1;1,1,2)={\displaystyle \frac{c_{1}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}{1-c_{2}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}}$$

(37)

for $c_1>0,c_2<0$,

$$V_{0,d_1,d_2}^{\left(1\right)}({\xi }_2;1,1,2)={\displaystyle \frac{d_{1}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}{1-d_{2}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}}$$

(38)

for $d_1>0,d_2<0$

$$V_{0,c_1,c_2}^{\left(2\right)}({\xi }_1;1,1,2)=-{\displaystyle \frac{c_{1}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}{1+c_{2}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}}$$

(39)

for $c_1<0,c_2>0$, and

$$V_{0,d_1,d_2}^{\left(2\right)}({\xi }_2;1,1,2)=-{\displaystyle \frac{d_{1}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}{1+d_{2}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}}$$

(40)

for $d_1<0,d_2>0$. The particular variants of ODEs of Bernoulli are

$${\displaystyle \frac{d{f}_1}{d{\xi }_1}}=c_1{f}_1+c_2{f}_1^2,{\displaystyle \frac{d{f}_2}{d{\xi }_2}}=d_1{f}_2+d_2{f}_2^2$$

(41)

Lemma 4. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (41). The application of the extended SEsM with simple equations of kind (41) reduces Equation (1) to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}\hfill & & -{a_1}^3+2a_1{c_2}^2=0\hfill \\ & & {\omega }_1{a}_1{c}_2+\beta {a_1}^2+3a_1{c}_1{c}_2+{a_1}^2-3a_0{a_1}^2=0\hfill \\ & & -a_1\alpha {b_1}^2+\sigma a_1\alpha {b_2}^2=0\hfill \\ & & -a_1{b}_1-\sigma a_1{b}_1-2\sigma a_1\alpha b_0{b}_1-2a_1\alpha b_0{b}_1=0\hfill \\ & & 2\beta a_0{a}_1-\sigma a_1\alpha {b_0}^2-\sigma a_1{b}_0-a_1{b}_0-\beta a_1\hfill \\ & & +2a_0{a}_1+{\omega }_1{a}_1{c}_1-a_1\alpha {b_0}^2+a_1{c_1}^2-3{a_0}^2{a}_1=0\hfill \\ & & {\mu }_3{b_1}^3-2\kappa b_1{d_2}^2=0\hfill \\ & & -a_0\alpha {b_1}^2+3{\mu }_3{b}_0{b_1}^2-\sigma a_0\alpha {b_1}^2-{\omega }_2{b}_1{d}_2-3\kappa b_1{d}_1{d}_2+{\mu }_2{b_1}^2=0\hfill \\ & & -{\omega }_2{b}_1{d}_1-a_0{b}_1+{\mu }_1{b}_1-2a_0\alpha b_0{b}_1-\kappa b_1{d_1}^2+2{\mu }_2{b}_0{b}_1\hfill \\ & & +3{\mu }_3{b_0}^2{b}_1-\sigma a_0{b}_1-2\sigma a_0\alpha b_0{b}_1=0\hfill \\ & & -\beta a_0+{a_0}^2+\beta {a_0}^2-a_0{b}_0+{\mu }_1{b}_0+{\mu }_2{b_0}^2+{\mu }_3{b_0}^3\hfill \\ & & -\sigma a_0{b}_0-a_0\alpha {b_0}^2-\sigma a_0\alpha {b_0}^2-{a_0}^3=0\hfill \end{array}$$

(42)

Proof. 

The application of the extended SEsM to Equation (1) means that we have to substitute Equation (7) for the case $n_1=n_2=1$, together with Equations (8) and (41) into Equation (1). Then, we have to equate to 0 the coefficients of the obtained polynomials of $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ in Equation (1). This leads us to the algebraic system presented in Equation (42). □

The combination of different solutions of the ODEs of Bernoulli, depending on the values of their coefficients, leads to obtaining several types of analytical solutions of Equation (1).

Proposition 4. 

The exact solutions of Equation (1) are

$$\begin{array}{ccc}\hfill & & u^{\left(7\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\alpha +{\Omega }_6}{\alpha }}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}{c}_1+\sqrt{2}{c}_2{V}_{0,c_1,c_2}^{\left(1\right)}({\xi }_1;1,1,2),\hfill \\ & & v^{\left(7\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_2}{{\mu }_3}}{V}_{0,d_1,d_2}^{\left(1\right)}({\xi }_2;1,1,2),\hfill \end{array}$$

(43)

where

$$V_{0,c_1,c_2}^{\left(1\right)}({\xi }_1;1,1,2)={\displaystyle \frac{c_{1}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}{1-c_{2}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}}$$

for $c_1>0,c_2<0$,

$$V_{0,d_1,d_2}^{\left(1\right)}({\xi }_2;1,1,2)={\displaystyle \frac{d_{1}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}{1-d_{2}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}}$$

for $d_1>0,d_2<0$, with ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}}$ and ${\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

$$\begin{array}{ccc}\hfill & & u^{\left(8\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\alpha +{\Omega }_6}{\alpha }}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}{c}_1+\sqrt{2}{c}_2{V}_{0,c_1,c_2}^{\left(2\right)}({\xi }_1;1,1,2),\hfill \\ & & v^{\left(8\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_2}{{\mu }_3}}{V}_{0,d_1,d_2}^{\left(2\right)}({\xi }_2;1,1,2),\hfill \end{array}$$

(44)

where

$$V_{0,c_1,c_2}^{\left(2\right)}({\xi }_1;1,1,2)=-{\displaystyle \frac{c_{1}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}{1+c_{2}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}}$$

for $c_1<0,c_2>0$,

$$V_{0,d_1,d_2}^{\left(2\right)}({\xi }_2;1,1,2)=-{\displaystyle \frac{d_{1}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}{1+d_{2}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}}$$

for $d_1<0,d_2>0$, with ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}}$ and ${\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

$$\begin{array}{ccc}\hfill & & u^{\left(9\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\alpha +{\Omega }_6}{\alpha }}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}{c}_1+\sqrt{2}{c}_2{V}_{0,c_1,c_2}^{\left(1\right)}({\xi }_1;1,1,2),\hfill \\ & & v^{\left(9\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_2}{{\mu }_3}}{V}_{0,d_1,d_2}^{\left(2\right)}({\xi }_2;1,1,2),\hfill \end{array}$$

(45)

where

$$V_{0,c_1,c_2}^{\left(1\right)}({\xi }_1;1,1,2)={\displaystyle \frac{c_{1}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}{1-c_{2}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}}$$

for $c_1>0,c_2<0$,

$$V_{0,d_1,d_2}^{\left(2\right)}({\xi }_2;1,1,2)=-{\displaystyle \frac{d_{1}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}{1+d_{2}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}}$$

for $d_1<0,d_2>0$, with ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}}$ and ${\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

$$\begin{array}{ccc}\hfill & & u^{\left(10\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\alpha +{\Omega }_6}{\alpha }}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}{c}_1+\sqrt{2}{c}_2{V}_{0,c_1,c_2}^{\left(2\right)}({\xi }_1;1,1,2),\hfill \\ & & v^{\left(10\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_1\kappa }{d}_2}{{\mu }_3}}{V}_{0,d_1,d_2}^{\left(1\right)}({\xi }_2;1,1,2),\hfill \end{array}$$

(46)

where

$$V_{0,c_1,c_2}^{\left(2\right)}({\xi }_1;1,1,2)=-{\displaystyle \frac{c_{1}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}{1+c_{2}exp\left[c_1({\xi }_1+{\xi }_{01})\right]}}$$

for $c_1<0,c_2>0$,

$$V_{0,d_1,d_2}^{\left(1\right)}({\xi }_2;1,1,2)={\displaystyle \frac{d_{1}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}{1-d_{2}exp\left[d_1({\xi }_2+{\xi }_{02})\right]}}$$

for $d_1>0,d_2<0$ and ${\xi }_1=x-{\displaystyle \frac{{\omega }_1{t}^{\varphi }}{\Gamma (1+\varphi )}},{\xi }_2=x-{\displaystyle \frac{{\omega }_2{t}^{\varphi }}{\Gamma (1+\varphi )}}$.

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of Equation (1) are of kind (41). We apply the extended SEsM to Equation (1) on the basis of Equation (7) for the case $n_1=n_2=1$, and Equations (8) and (41). As a result, we obtain the system of nonlinear algebraic equations presented in Equation (42). A non–trivial solution of this system is

$$\begin{array}{c}\hfill a_0={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\alpha +{\Omega }_6}{\alpha }}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}{c}_1,a_1=\sqrt{2}{c}_2,b_0=-{\displaystyle \frac{1}{2\alpha }},b_1=-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_2}{{\mu }_3}}\\ \hfill \beta ={\displaystyle \frac{1}{2}}{\displaystyle \frac{\alpha +{\Omega }_6}{\alpha }},\sigma =-1,{\mu }_2={\displaystyle \frac{1}{4}}{\displaystyle \frac{{\mu }_3\left(-3{\mu }_3+6\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_1\alpha -8\kappa {d_1}^2{\alpha }^2-4{\mu }_1{\alpha }^2\right)}{\left(\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_1\alpha -{\mu }_3\right)\alpha }}\\ \hfill {\omega }_2=-{\displaystyle \frac{1}{4}}{\displaystyle \frac{3\sqrt{2}\sqrt{{\mu }_3\kappa }{\mu }_3-12{\mu }_3\kappa \alpha d_1+4\kappa \sqrt{2}\sqrt{{\mu }_3\kappa }{d_1}^2{\alpha }^2-4{\mu }_1\sqrt{2}\sqrt{{\mu }_3\kappa }{\alpha }^2}{\left(\sqrt{2}\sqrt{{\mu }_3\kappa }{d}_1\alpha -{\mu }_3\right)\alpha }}\end{array}$$

(47)

where

$${\Omega }_6=\sqrt{-3{\alpha }^2+2{{\omega }_1}^2{\alpha }^2+6{c_1}^2{\alpha }^2}$$

Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (47) into (36) together with (37)–(40) leads to the solutions (43)–(46) of (1). □

For the special case $c_1=-1,c_2=1$ and $d_1=-1,d_2=1$, the solution (36) reduces to

$$\begin{array}{c}\hfill u\left({\xi }_1\right)=a_0+a_1{V}_{0,-1,1}({\xi }_1;1,1,2),\\ \hfill v\left({\xi }_2\right)=b_0+b_1{V}_{0,-1,1}({\xi }_2;1,1,2),\end{array}$$

(48)

where the special functions $V_{0,-1,1}({\xi }_1;1,1,2)$ and $V_{0,-1,1}({\xi }_2;1,1,2)$ are the solutions of particular variants of an ODE of Riccati [39].

$$V_{0,-1,1}({\xi }_1;1,1,2)={\displaystyle \frac{1}{1+exp({\xi }_1+{\xi }_{01})}},V_{0,-1,1}({\xi }_2;1,1,2)={\displaystyle \frac{1}{1+exp({\xi }_2+{\xi }_{02})}}$$

(49)

The particular variants of an ODEs of Riccati are

$${\displaystyle \frac{d{f}_1}{d{\xi }_1}}=f_1^2-f_1,{\displaystyle \frac{d{f}_2}{d{\xi }_2}}=f_2^2-f_2$$

(50)

Lemma 5. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (50). The application of the extended SEsM with simple equations of kind (50) reduces Equation (1) to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}\hfill & & 2a_1-{a_1}^3=0\hfill \\ & & -3a_1+{\omega }_1{a}_1-3a_0{a_1}^2+\beta {a_1}^2+{a_1}^2=0\hfill \\ & & -\sigma a_1\alpha {b_1}^2+a_1\alpha {b_2}^2=0\hfill \\ & & -a_1{b}_1-2a_1\alpha b_0{b}_1-2\sigma a_1\alpha b_0{b}_1-\sigma a_1{b}_1=0\hfill \\ & & -\sigma a_1{b}_0-\sigma a_1\alpha {b_0}^2+a_1-a_1{b}_0+2a_0{a}_1-{\omega }_1{a}_1-a_1\alpha {b_0}^2-\beta a_1+2\beta a_0{a}_1-3{a_0}^2{a}_1=0\hfill \\ & & {\mu }_3{b_1}^3-2\kappa b_1=0\hfill \\ & & -{\omega }_2{b}_1+{\mu }_2{b_1}^2-\sigma a_0\alpha {b_1}^2-a_0\alpha {b_1}^2+3{\mu }_3{b}_0{b_1}^2+3\kappa b_1=0\hfill \\ & & -2a_0\alpha b_0{b}_1+{\omega }_2{b}_1+2{\mu }_2{b}_0{b}_1-a_0{b}_1-2\sigma a_0\alpha b_0{b}_1+3{\mu }_3{b_0}^2{b}_1-\sigma a_0{b}_1+{\mu }_1{b}_1-\kappa b_1=0\hfill \\ & & -a_0{b}_0+{\mu }_1{b}_0+{\mu }_2{b_0}^2+{\mu }_3{b_0}^3+\beta {a_0}^2-\sigma a_0{b}_0-\beta a_0-{a_0}^3-a_0\alpha {b_0}^2+{a_0}^2-\sigma a_0\alpha {b_0}^2=0\hfill \end{array}$$

(51)

Proof. 

The application of the extended SEsM to Equation (1) means that we have to substitute (7) for the case $n_1=n_2=1$, together with (8) and (50) into (1). Then, we have to equate to 0 the coefficients of the obtained polynomials of $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ in (1). This leads us to the algebraic system presented in Equation (51). □

Proposition 5. 

An exact solution of Equation (1) is

$$\begin{array}{ccc}\hfill & & u^{\left(11\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{3}}-{\displaystyle \frac{1}{6}}\sqrt{-2\beta +2{\beta }^2-1}\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}+{\displaystyle \frac{1}{3}}\beta +\sqrt{2}{V}_{0,-1,1}({\xi }_1;1,1,2)\hfill \\ & & v^{\left(11\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}+{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{0,-1,1}({\xi }_2;1,1,2)\hfill \end{array}$$

(52)

where

$$V_{0,-1,1}({\xi }_1;1,1,2)={\displaystyle \frac{1}{1+exp({\xi }_1+{\xi }_{01})}},V_{0,-1,1}({\xi }_2;1,1,2)={\displaystyle \frac{1}{1+exp({\xi }_2+{\xi }_{02})}}$$

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

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of Equation (1) are of kind (50). We apply the extended SEsM to Equation (1) on the basis of Equation (7) for the case $n_1=n_2=1$, and Equations (8) and (50). As a result, we obtain the system of nonlinear algebraic equations presented in Equation (51). A non-trivial solution of this system is

$$\begin{array}{ccc}\hfill & & a_0={\displaystyle \frac{1}{3}}-{\displaystyle \frac{1}{6}}\sqrt{-2\beta +2{\beta }^2-1}\sqrt{2}+{\displaystyle \frac{1}{2}}\sqrt{2}+{\displaystyle \frac{1}{3}}\beta ,a_1=\sqrt{2},b_0=-{\displaystyle \frac{1}{2\alpha }},b_1=+{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}},\hfill \\ & & {\omega }_1=-\sqrt{-2\beta +2{\beta }^2-1},\sigma ={\displaystyle \frac{2}{3}}{\omega }_1\alpha -{\displaystyle \frac{4}{3}}{\beta }^2\alpha +{\displaystyle \frac{2}{3}}\alpha -1+{\displaystyle \frac{4}{3}}\beta \alpha ,\hfill \\ & & {\omega }_2={\displaystyle \frac{1}{54}}{\displaystyle \frac{{\mu }_3}{\left(2\sqrt{2}\sqrt{{\mu }_3\kappa }\alpha -{\mu }_3\right)\alpha \left(\sqrt{2}\sqrt{{\mu }_3\kappa }\alpha -{\mu }_3\right)}}(-270\kappa \sqrt{2}\sqrt{{\mu }_3\kappa }{\alpha }^2+216{\mu }_3\kappa \alpha \hfill \\ & & -27\sqrt{2}\sqrt{{\mu }_3\kappa }{\mu }_3+48\kappa \beta {\alpha }^4-24\sqrt{2}\sqrt{{\mu }_3\kappa }\beta {\alpha }^3-80\kappa \sqrt{2}{\omega }_1{\alpha }^4+80\sqrt{{\mu }_3\kappa }{\omega }_1{\alpha }^3-32\kappa {\alpha }^4\hfill \\ & & +16\sqrt{2}\sqrt{{\mu }_3\kappa }{\alpha }^3+16\kappa \sqrt{2}{\omega }_1{\beta }^2{\alpha }^4-16{\omega }_1{\beta }^2{\alpha }^3-32\kappa {\beta }^3{\alpha }^4+16\sqrt{2}\sqrt{{\mu }_3\kappa }{\beta }^3{\alpha }^3\hfill \\ & & -16\kappa \sqrt{2}{\omega }_1\beta {\alpha }^4+16\sqrt{{\mu }_3\kappa }{\omega }_1\beta {\alpha }^3+48\kappa {\beta }^2{\alpha }^4-24\sqrt{2}\sqrt{{\mu }_3\kappa }{\beta }^2{\alpha }^3+216{\alpha }^3{\kappa }^2)\hfill \\ & & {\mu }_1={\displaystyle \frac{1}{108}}{\displaystyle \frac{{\Omega }_7}{{\alpha }^2\left(2\sqrt{2}\sqrt{{\mu }_3\kappa }\alpha -{\mu }_3\right)}},{\mu }_2=-{\displaystyle \frac{1}{27}}{\displaystyle \frac{{\mu }_3{\Omega }_7}{\left(2\sqrt{2}\sqrt{{\mu }_3\kappa }\alpha -{\mu }_3\right)\left(\sqrt{2}\sqrt{{\mu }_3\kappa }\alpha -{\mu }_3\right)\alpha }}\hfill \end{array}$$

(53)

where

$${\Omega }_7=108\sqrt{2}\sqrt{{\mu }_3\kappa }{\mu }_3\alpha -27{{\mu }_3}^2-32{\mu }_3{\beta }^3{\alpha }^3-48{\beta }^2{\alpha }^4\sqrt{2}\sqrt{{\mu }_3\kappa }+48{\mu }_3{\beta }^2{\alpha }^3$$

Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (53) into (48), together with (49), leads to solution (52) of (1). □

3.2.3. Variant 3: When $c_0\ne 0,c_1=0,c_2\ne 0$ and $d_0\ne 0,d_1=0,d_2\ne 0$ in Equation (29)

For this case, the general solution of Equation (1) can be presented as

$$\begin{array}{c}\hfill u\left({\xi }_1\right)=a_0+a_1{V}_{c_0,0,c_2}({\xi }_1;1,1,2),\\ \hfill v\left({\xi }_2\right)=b_0+b_1{V}_{d_0,0,d_2}({\xi }_2;1,1,2),\end{array}$$

(54)

where $V_{c_0,0,c_2}({\xi }_1;1,1,2)$ and $V_{d_0,0,d_2}({\xi }_2;1,1,2)$ are particular solutions of the reduced Riccati equations, known also as $tanh$-function equations [39].

$$\begin{array}{c}\hfill V_{c_0,0,c_2}({\xi }_1;1,1,2)={\displaystyle \frac{c_0}{c_2}}tanh\left[c_2{c}_0({\xi }_1+{\xi }_{01})\right],\\ \hfill V_{d_0,0,d_2}({\xi }_2;1,1,2)={\displaystyle \frac{d_0}{d_2}}tanh\left[d_2{d}_0({\xi }_2+{\xi }_{02})\right],\end{array}$$

(55)

The $tanh$–function equations are

$${\displaystyle \frac{d{f}_1}{d{\xi }_1}}=c_0+c_2{f}_1^2,{\displaystyle \frac{d{f}_2}{d{\xi }_2}}=d_0+d_2{f}_2^2$$

(56)

Lemma 6. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (56). The application of the extended SEsM with simple equations of kind (56) reduces Equation (1) to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}\hfill & & -{a_1}^3+2a_1{c_2}^2=0\hfill \\ & & \beta {a_1}^2+{a_1}^2+{\omega }_1{a}_1{c}_2-3a_0{a_1}^2=0\hfill \\ & & -\sigma a_1\alpha {b_1}^2+a_1\alpha {b_2}^2=0\hfill \\ & & -a_1{b}_1-\sigma a_1{b}_1-2\sigma a_1\alpha b_0{b}_1-2a_1\alpha b_0{b}_1=0\hfill \\ & & -\sigma a_1{b}_0-\sigma a_1\alpha {b_0}^2+2a_0{a}_1-a_1\alpha {b_0}^2-\beta a_1\hfill \\ & & +2\beta a_0{a}_1+2a_1{c}_2{c}_0-3{a_0}^2{a}_1-a_1{b}_0=0\hfill \end{array}$$

(57)

$$\begin{array}{ccc}& & {\mu }_3{b_1}^3-2\kappa b_1{d_2}^2=0\hfill \\ & & -\sigma a_0\alpha {b_1}^2+3{\mu }_3{b}_0{b_1}^2-a_0\alpha {b_1}^2-{\omega }_2{b}_1{d}_2+{\mu }_2{b_1}^2=0\hfill \\ & & -2\kappa b_1{d}_2{d}_0-\sigma a_0{b}_1+2{\mu }_2{b}_0{b}_1+3{\mu }_3{b_0}^2{b}_1\hfill \\ & & +{\mu }_1{b}_1-a_0{b}_1-2a_0\alpha b_0{b}_1-2\sigma a_0\alpha b_0{b}_1=0\hfill \\ & & {\omega }_1{a}_1{c}_0-\beta a_0+{a_0}^2+\beta {a_0}^2-a_0{b}_0-\sigma a_0\alpha {b_0}^2+{\mu }_2{b_0}^2\hfill \\ & & +{\mu }_3{b_0}^3-a_0\alpha {b_0}^2-\sigma a_0{b}_0-{a_0}^3-{\omega }_2{b}_1{d}_0+{\mu }_1{b}_0=0\hfill \end{array}$$

Proof. 

The application of the extended SEsM to Equation (1) means that we have to substitute Equation (7) for the case $n_1=n_2=1$, together with Equations (8) and (56) into Equation (1). Then, we have to equate to 0 the coefficients of the obtained polynomials of $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ in Equation (1). This leads us to the algebraic system given in (57). □

Proposition 6. 

An exact solution of Equation (1) is

$$\begin{array}{ccc}\hfill & & u^{\left(12\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\alpha +{\Omega }_7}{\alpha }}+\sqrt{2}{c}_2{V}_{c_0,0,c_2}({\xi }_1;1,1,2),\hfill \\ & & v^{\left(12\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}+{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_1\kappa }{d}_2}{{\mu }_3}}{V}_{d_0,0,d_2}({\xi }_2;1,1,2),\hfill \end{array}$$

(58)

where

$$\begin{array}{c}\hfill V_{c_0,0,c_2}({\xi }_1;1,1,2)={\displaystyle \frac{c_0}{c_2}}tanh\left[c_2{c}_0({\xi }_1+{\xi }_{01})\right],\\ \hfill V_{d_0,0,d_2}({\xi }_2;1,1,2)={\displaystyle \frac{d_0}{d_2}}tanh\left[d_2{d}_0({\xi }_2+{\xi }_{02})\right],\end{array}$$

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

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of Equation (1) are of kind (56). We apply the extended SEsM to Equation (1) on the basis of Equation (7) for the case $n_1=n_2=1$, and Equations (8) and (56). As a result, we obtain the system of nonlinear algebraic Equation (57). A non-trivial solution of this system is

$$\begin{array}{ccc}\hfill & & a_0={\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}{\omega }_1\sqrt{2}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{\alpha +{\Omega }_8}{\alpha }},a_1=\sqrt{2}{c}_2,b_0=-{\displaystyle \frac{1}{2\alpha }},b_1={\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_1\kappa }{d}_2}{{\mu }_3}},\beta ={\displaystyle \frac{1}{2}}{\displaystyle \frac{\alpha +{\Omega }_8}{\alpha }},\hfill \\ & & {\mu }_2=-{\displaystyle \frac{1}{18}}{\displaystyle \frac{9{\mu }_1{\alpha }^2-18{\mu }_3+\sqrt{3}\alpha \sqrt{-12{\Omega }_8{\mu }_3-27{\alpha }^2{{\mu }_1}^2+72\alpha {\omega }_1\sqrt{2}{c}_2{c}_0{\mu }_3-24{\Omega }_8{c}_2{c}_0{\mu }_3}}{\alpha }}\hfill \\ & & {\omega }_2={\displaystyle \frac{1}{4}}{\mu }_1{\alpha }^2-{\displaystyle \frac{1}{8}}{\mu }_3+{\displaystyle \frac{1}{36}}{\displaystyle \frac{\left(\sqrt{36\alpha {\Omega }_7{\mu }_3+81{\mu }_1^2{\alpha }^2-6{{\omega }_1}^3\sqrt{2}\alpha {\mu }_3+72{\Omega }_8{c}_2{c}_0{\mu }_3}\right)}{{\mu }_3}},\hfill \end{array}$$

(59)

where

$${\Omega }_8=\sqrt{-3{\alpha }^2-24c_2{c}_0{\alpha }^2+2{{\omega }_1}^2{\alpha }^2}$$

Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (59) into (54), together with (55), leads to the solution (58) of (1). □

4. Exact Solutions of Time-Fractional Diffusive Predator–Prey System Incorporating Allee Effect Using Standard Traveling Wave Transformation

For this scenario, System (1) is defined in sense of the modified Riemann–Liouville derivatives. The general solution of Equation (1) is presented by Equation (7). We introduce the following traveling wave transformations:

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

(60)

The simple equations used are as follows [32,33,34]:

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

(61)

The balance equations derived are the same as those shown in Section 3. Thus, we express Equation (7) for the case $n_1=n_2=1$ by special functions, as follows:

$$u\left({\xi }_1\right)=a_0+a_1{V}_{c_0,0,1}({\xi }_1,\varphi ;2),v\left({\xi }_2\right)=b_0+b_1{V}_{d_0,0,1}({\xi }_2,\varphi ;2),$$

(62)

where $V_{c_0,0,1}({\xi }_1,\varphi ;2)$ and $V_{d_0,0,1}({\xi }_2,\varphi ;2)$ are appropriate solutions of Equation (61) taken from [32,33,34]:

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

(63)

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

(64)

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

(65)

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

(66)

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

(67)

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

(68)

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

(69)

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

(70)

The generalized trigonometric and hyperbolic functions presented in Equations (63)–(70) are defined as follows:

$$\begin{array}{ccc}\hfill & & 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 \end{array}$$

(71)

$$\begin{array}{ccc}& & 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}$$

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.

Lemma 7. 

Let us consider the case where the simple equations used for finding exact solutions of (1) are of kind (61). The application of the extended SEsM with simple equations of kind (61) reduces Equation (1) to the following system of nonlinear algebraic equations:

$$\begin{array}{ccc}\hfill & & -{a_1}^3+2a_1=0\hfill \\ & & \beta {a_1}^2+{a_1}^2+{\omega }_1^{\varphi }{a}_1-3a_0{a_1}^2=0\hfill \\ & & -\sigma a_1\alpha {b_1}^2+a_1\alpha {b_2}^2=0\hfill \\ & & -a_1{b}_1-\sigma a_1{b}_1-2\sigma a_1\alpha b_0{b}_1-2a_1\alpha b_0{b}_1=0\hfill \\ & & -\sigma a_1{b}_0-\sigma a_1\alpha {b_0}^2+2a_0{a}_1-a_1\alpha {b_0}^2-\beta a_1\hfill \\ & & +2\beta a_0{a}_1+2a_1{c}_0-3{a_0}^2{a}_1-a_1{b}_0=0\hfill \\ & & {\mu }_3{b_1}^3-2\kappa b_1=0\hfill \\ & & -\sigma a_0\alpha {b_1}^2+3{\mu }_3{b}_0{b_1}^2-a_0\alpha {b_1}^2-{\omega }_2^{\varphi }{b}_1+{\mu }_2{b_1}^2=0\hfill \\ & & -2\kappa b_1{d}_0-\sigma a_0{b}_1+2{\mu }_2{b}_0{b}_1+3{\mu }_3{b_0}^2{b}_1\hfill \\ & & +{\mu }_1{b}_1-a_0{b}_1-2a_0\alpha b_0{b}_1-2\sigma a_0\alpha b_0{b}_1=0\hfill \\ & & {\omega }_1^{\varphi }{a}_1{c}_0-\beta a_0+{a_0}^2+\beta {a_0}^2-a_0{b}_0-\sigma a_0\alpha {b_0}^2+{\mu }_2{b_0}^2\hfill \\ & & +{\mu }_3{b_0}^3-a_0\alpha {b_0}^2-\sigma a_0{b}_0-{a_0}^3-{\omega }_2^{\varphi }{b}_1{d}_0+{\mu }_1{b}_0=0\hfill \end{array}$$

(72)

Proof. 

The application of the extended SEsM to Equation (1) means that we have to substitute Equation (7) for the case $n_1=n_2=1$, s (60) and (61) into Equation (1). Then, we have to equate to 0 the coefficients of the obtained polynomials of $f_1\left({\xi }_1\right)$ and $f_2\left({\xi }_2\right)$ in Equation (1). This leads us to the algebraic system (72). □

The combination of different solutions of the fractional ODEs of Riccati, depending on the numerical values of their coefficients, leads to a wide range of different analytical solutions of Equation (1). Below, all possible exact solutions of System (1) are presented.

4.1. Case 1: When $c_0<0$ and $d_0<0$ in Equation (61)

Let us consider the case where $c_0<0$ and $d_0<0$ in Equation (61).

Proposition 7. 

The exact solutions of Equation (1) are

$$\begin{array}{ccc}\hfill & & u^{\left(13\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(13\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(73)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}tanh[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}tanh[\sqrt{-d_0}{\xi }_2,\varphi ],\end{array}$$

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

$$\begin{array}{ccc}\hfill & & u^{\left(14\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(14\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(74)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}coth[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;1,2)=-\sqrt{-d_0}coth[\sqrt{-d_0}{\xi }_2,\varphi ],\end{array}$$

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

$$\begin{array}{ccc}\hfill & & u^{\left(15\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(15\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(75)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}tanh[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}coth[\sqrt{-d_0}{\xi }_2,\varphi ],\end{array}$$

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

$$\begin{array}{ccc}\hfill & & u^{\left(16\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(16\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(76)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}coth[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}tanh[\sqrt{-d_0}{\xi }_2,\varphi ],\end{array}$$

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

Proof. 

Let us consider the case when the simple equations used for finding exact solutions of Equation (1) are of kind (61), where $c_0<0$ and $d_0<0$. We apply the extended SEsM to Equation (1) on the basis of Equation (7) in the case $n_1=n_2=1$, and Equations (60) and (61). As a result, we obtain the system of nonlinear algebraic Equation (72). A non-trivial solution of this system is

$$\begin{array}{ccc}\hfill & & a_0={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2},a_1=-\sqrt{2},b_0=-{\displaystyle \frac{1}{2\alpha }},b_1=-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}},\hfill \\ & & \beta ={\displaystyle \frac{1}{2}}(1+{\Omega }_9),\sigma =-1,{\mu }_2=-{\displaystyle \frac{{\mu }_3}{27{\Omega }_{10}}}(-27{\Omega }_9+{\Omega }_9{\Omega }_{11}+{\Omega }_{12}),\hfill \\ & & {\mu }_1=-{\displaystyle \frac{1}{\alpha {\Omega }_{10}}}(-432{\Omega }_{10}+4{\Omega }_9{\Omega }_{11}+4{\Omega }_{12}),{\omega }_2^{\varphi }={\displaystyle \frac{\sqrt{{\mu }_3\kappa }}{54{\Omega }_{10}}}(27{\Omega }_{10}-2{\Omega }_9{\Omega }_{11}+2{\Omega }_{12}),\hfill \end{array}$$

(77)

where

$$\begin{array}{ccc}& & {\Omega }_9=\sqrt{-3-24c_0+2{{\omega }_1}^{2\varphi }},{\Omega }_{10}={\mu }_3+8\kappa d_0{\alpha }^2,\hfill \\ & & {\Omega }_{11}=24{\alpha }^3{c}_0-2{\alpha }^3{{\omega }_1^{\varphi }}^2+12{\alpha }^3,{\Omega }_{12}=2{{\omega }_1^{\varphi }}^3\sqrt{2}{\alpha }^3+72{\omega }_1^{\varphi }\sqrt{2}{c}_0{\alpha }^3\hfill \end{array}$$

Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (77) into (62), together with (63)–(66), leads to the solutions (73)–(76) of (1). □

4.2. Case 2: When $c_0>0$ and $d_0>0$ in Equation (61)

Let us consider the case where $c_0>0$ and $d_0>0$ in Equation (61).

Proposition 8. 

The exact solutions of Equation (1) are

$$\begin{array}{ccc}\hfill & & u^{\left(17\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(17\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(78)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)=\sqrt{c_0}tan[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)=\sqrt{d_0}tan[\sqrt{d_0}{\xi }_2,\varphi ],\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne {\displaystyle \frac{\pi }{2}}+n\pi$ and ${\xi }_2=x-{\omega }_{2}t\ne {\displaystyle \frac{\pi }{2}}+n\pi ,\pi \in Z$;

$$\begin{array}{ccc}\hfill & & u^{\left(18\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}^{{\left(4\right)}_{c_0,0,1}({\xi }_1,\varphi ;2)},\hfill \\ & & v^{\left(18\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}^{{\left(4\right)}_{d_0,0,1}({\xi }_2,\varphi ;2)}\hfill \end{array}$$

(79)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}cot[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}cot[\sqrt{d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne n\pi$ and ${\xi }_2=x-{\omega }_{2}t\ne n\pi ,\pi \in Z$;

$$\begin{array}{ccc}\hfill & & u^{\left(19\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(19\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(80)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)=\sqrt{c_0}tan[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}cot[\sqrt{d_0}{\xi }_2,\varphi ],\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne {\displaystyle \frac{\pi }{2}}+n\pi$ and ${\xi }_2=x-{\omega }_{2}t\ne n\pi ,\pi \in Z$;

$$\begin{array}{ccc}\hfill & & u^{\left(20\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(20\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(81)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}cot[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)=\sqrt{d_0}tan[\sqrt{d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne n\pi$ and ${\xi }_2=x-{\omega }_{2}t\ne {\displaystyle \frac{\pi }{2}}+n\pi ,\pi \in Z$.

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of Equation (1) are of kind (61), where $c_0>0$ and $d_0>0$. We apply the extended SEsM to Equation (1) on the basis of Equation (7) for the case $n_1=n_2=1$, and Equations (60) and (61). As a result, we obtain the system of nonlinear algebraic Equation (72). A non-trivial solution of this system is (77). Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (77) into (62), together with (67)–(70), leads to the solutions (78)–(81) of (1). □

4.3. Case 3: When $c_0<0$ and $d_0>0$ in Equations (61)

Let us consider the case where $c_0<0$ and $d_0>0$ in Equation (61).

Proposition 9. 

The exact solutions of Equation (1) are

$$\begin{array}{ccc}\hfill & & u^{\left(21\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(22\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}^{{\left(3\right)}_{d_0,0,1}({\xi }_2,\varphi ;2)}\hfill \end{array}$$

(82)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}tanh[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)=\sqrt{d_0}tan[\sqrt{d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t$ and ${\xi }_2=x-{\omega }_{2}t\ne {\displaystyle \frac{\pi }{2}}+n\pi ,\pi \in Z$;

$$\begin{array}{ccc}\hfill & & u^{\left(22\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(22\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(83)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(1\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}tanh[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}cot[\sqrt{d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t$ and ${\xi }_2=x-{\omega }_{2}t\ne n\pi ,\pi \in Z$;

$$\begin{array}{ccc}\hfill & & u^{\left(23\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(23\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(84)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}coth[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(3\right)}({\xi }_2,\varphi ;2)=\sqrt{d_0}tan[\sqrt{d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne 0$ and ${\xi }_2=x-{\omega }_{2}t\ne {\displaystyle \frac{\pi }{2}}+n\pi ,\pi \in Z$;

$$\begin{array}{ccc}\hfill & & u^{\left(24\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(24\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(85)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(2\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}coth[\sqrt{-c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(4\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}cot[\sqrt{d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne 0$ and ${\xi }_2=x-{\omega }_{2}t\ne n\pi ,\pi \in Z$.

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of Equation (1) are of kind (61), where $c_0<0$ and $d_0>0$. We apply the extended SEsM to Equation (1) on the basis of Equation (7) for the case $n_1=n_2=1$, and Equations (60) and (61). As a result, we obtain the system of nonlinear algebraic Equation (72). A non-trivial solution of this system is (77). Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (77) into (62), together with Equations (63), (65), (68) and (70), leads to the solutions (82)–(85) of (1). □

4.4. Case 4: When $c_0>0$ and $d_0<0$ in Equation (61)

Let us consider the case where $c_0>0$ and $d_0<0$ in Equation (61).

Proposition 10. 

The exact solutions of Equation (1) are

$$\begin{array}{ccc}\hfill & & u^{\left(25\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(25\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(86)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)=\sqrt{c_0}tan[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}tanh[\sqrt{-d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne {\displaystyle \frac{\pi }{2}}+n\pi$ and ${\xi }_2=x-{\omega }_{2}t$;

$$\begin{array}{ccc}\hfill & & u^{\left(26\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(26\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(87)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(3\right)}({\xi }_1,\varphi ;2)=\sqrt{c_0}tan[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}coth[\sqrt{-d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne {\displaystyle \frac{\pi }{2}}+n\pi$ and ${\xi }_2=x-{\omega }_{2}t\ne 0$;

$$\begin{array}{ccc}\hfill & & u^{\left(27\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(27\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(88)

where

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}cot[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(1\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}tanh[\sqrt{-d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne n\pi$ and ${\xi }_2=x-{\omega }_{2}t$;

$$\begin{array}{ccc}\hfill & & u^{\left(28\right)}\left({\xi }_1\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6}}{\Omega }_9-{\displaystyle \frac{1}{6}}{{\omega }_1}^{\varphi }\sqrt{2}-\sqrt{2}{V}_{c_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2),\hfill \\ & & v^{\left(28\right)}\left({\xi }_2\right)=-{\displaystyle \frac{1}{2\alpha }}-{\displaystyle \frac{\sqrt{2}\sqrt{{\mu }_3\kappa }}{{\mu }_3}}{V}_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)\hfill \end{array}$$

(89)

$$\begin{array}{c}\hfill V_{c_0,0,1}^{\left(4\right)}({\xi }_1,\varphi ;2)=-\sqrt{-c_0}cot[\sqrt{c_0}{\xi }_1,\varphi ],V_{d_0,0,1}^{\left(2\right)}({\xi }_2,\varphi ;2)=-\sqrt{-d_0}coth[\sqrt{-d_0}{\xi }_2,\varphi ]\end{array}$$

with ${\xi }_1=x-{\omega }_{1}t\ne n\pi$ and ${\xi }_2=x-{\omega }_{2}t\ne 0$.

Proof. 

Let us consider the case where the simple equations used for finding exact solutions of Equation (1) are of kind (61), where $c_0>0$ and $d_0<0$. We apply the extended SEsM to Equation (1) on the basis of Equation (7) for the case $n_1=n_2=1$, and Equations (60) and (61). As a result, we obtain the system of nonlinear algebraic Equation (72). A non-trivial solution of this system is (77). Substitution of $a_0$ and $a_1$ and $b_0$ and $b_1$ from (77) into (62), together with Equations (64), (66), (67) and (69), leads to the solutions (86)–(89) of (1). □

All analytical calculations in Section 3 and Section 4 have been carried out using the mathematical computer software Maple 15 (Waterloo Maple Inc., https://www.maplesoft.com, accessed on 15 August 2024). Thanks to the Maple software, the solutions of all algebraic systems given in Section 3 and Section 4 were obtained based on a stepwise mechanism, starting from the simplest equations of the systems of algebraic equations and then moving towards solving more complicated equations. I note that, by means of the computational software, many variants of solutions of the algebraic systems in Section 3 and Section 4 can be obtained. However, my goal is to show how by means of the extended SEsM, we can generate many different exact solutions of Equation (1) with only one set of their coefficients (i.e., one solution of the corresponding algebraic systems) by varying only the solutions of the simple equations used for each of the cases presented in Section 3 and Section 4.

5. Numerical Results and Discussions

It is obvious that in the studied time-fractional predator–prey system, the corresponding population waves can demonstrate different dynamics than that of the classic predator–prey model. This is due to the added model parameters, which account for the strength of the Allee effect. Moreover, the studied predator–prey interaction is particularly influenced by the fractional order. To illustrate the influence of each parameter mentioned above on the dynamic wave behavior of the studied model, I have numerically simulated only one of the obtained exact solutions of (1) (Equation (34) in Section 3), using the computer software Maple. This was carried out in order to be able to compare the different types of waves that can appear under the same conditions by changing only the numerical values of key model parameters. In Figure 1(a*,a**), Equation (34) is plotted for the case $\varphi =1$, i.e., for the case of the standard (integer-order) version of (1). In this case, $\alpha =0.5$ (The parameter $\alpha$ is (the non-dimensional Allee threshold of the predator, leading to changing the non-dimensional Allee threshold of the prey $\beta$ in (1).) As can be seen, bistable waves moving from a lower stable state to a higher stable state are observed for both prey and predator populations. In Figure 1(b*,b**), the same solution is plotted when $\varphi =0.55$ and $\alpha =0.05$. As can be seen from Figure 1(b*,b**), the memory effect arising from the fractional derivative order slows down the rate of change in population densities, as the Allee effect in a time-fractional predator–prey system can lead to significantly more complex wave dynamics in both populations.This is manifested in the simultaneous occurrence of two distinct waves:a kink-type wave and a single soliton, observed in the predator and the prey population densities. We note that such complexity cannot be captured by the standard diffusion-based predator–prey model. Finally, although classical periodic waves are not typical for the studied population model, due to the memory effect and the stabilizing influence of the Allee thresholds, the system may display also damped or quasi-periodic oscillations that diminish over time (see Figure 1(c*,c**)). The last illustrations are obtained by numerical simulations of another solution of (1) (Equation (43) in Section 3). In this case, when the initial population densities of prey and predator exceed the Allee threshold, oscillations can occur, as both populations are large enough to interact meaningfully. The Allee effect, however, can also lead to damped oscillations, especially if the initial conditions are close to the Allee threshold. In these cases, oscillations may begin, but they can gradually reduce in amplitude until the system settles into a stable state. In addition, the fractional order induces a memory effect, which leads to slowing the oscillation frequency and creating smoother, longer-lasting oscillations, as can be seen in Figure 1(c*,c**). Numerical examples given here have illustrated only a part of the complex wave dynamics that can be observed in this specific predator–prey interaction. Building on the analytical results presented in Section 3 and Section 4, it is interesting to see what types of waves can arise from simulations of other exact solutions of the model, as well as how the interaction between key model parameters affects these dynamics. However, this will be a subject of further investigation.

6. Conclusions

In this study, an extended version of the SEsM was proposed for finding analytical solutions of systems of FNPDEs, whose system variables exhibit different wave dynamics. By the extended SEsM algorithm, which includes different simple equations with different independent variables in the system solutions, this work aims to provide a more flexible and comprehensive framework for a more realistic understanding of the complex wave dynamics of real-world processes modeled by FNPDEs. In addition, the introduction of several traveling wave transformations to the SEsM algorithm allows for reducing FNPDEs to several types ODEs, which leads to significantly increasing the spectrum of obtained analytical solutions to such systems.

As an example of the application of the proposed extended version of the SEsM, the time-fractional diffuse predator–prey system, including the Allee effect, was considered. It was assumed that the populations of the predator and the prey in the predator–prey interaction model considered moved at different wave speeds. From an applied point of view, this is a quite common situation in nature (for example, in animal populations). Based on the extended SEsM, various traveling wave type solutions of the studied system were obtained, one of which was numerically simulated at different values of the fractional order of the system equations and the Allee thresholds. These findings open the door to further questions about the wave dynamics in the time-fractional predator–prey model. For example, of particular interest is the question of what types of waves may arise from simulations of other exact solutions of the model, as well as how the interplay between key parameters influences these dynamics. By conducting additional numerical simulations across a broader range of parameter values and solutions, a more comprehensive comparative analysis of the observed wave patterns can be achieved. Such an investigation, however, extends beyond the current scope and will be the focus of future in-depth research.

Furthermore, for more realistic modeling of such ecological systems, our future research plans will be aimed at extending current models by integrating age-structured dynamics into them, or incorporating more realistic spatial movement models, in line with recent studies [62,63,64,65]. On the other hand, extending the used methodology to multidimensional systems or using more complex fractional operators could also help to improve both the realism of modeled ecological systems and their future applicability.