Simple Equations Method (SEsM): An Effective Algorithm for Obtaining Exact Solutions of Nonlinear Differential Equations

Exact solutions of nonlinear differential equations are of great importance to the theory and practice of complex systems. The main point of this review article is to discuss a specific methodology for obtaining such exact solutions. The methodology is called the SEsM, or the Simple Equations Method. The article begins with a short overview of the literature connected to the methodology for obtaining exact solutions of nonlinear differential equations. This overview includes research on nonlinear waves, research on the methodology of the Inverse Scattering Transform method, and the method of Hirota, as well as some of the nonlinear equations studied by these methods. The overview continues with articles devoted to the phenomena described by the exact solutions of the nonlinear differential equations and articles about mathematical results connected to the methodology for obtaining such exact solutions. Several articles devoted to the numerical study of nonlinear waves are mentioned. Then, the approach to the SEsM is described starting from the Hopf–Cole transformation, the research of Kudryashov on the Method of the Simplest Equation, the approach to the Modified Method of the Simplest Equation, and the development of this methodology towards the SEsM. The description of the algorithm of the SEsM begins with the transformations that convert the nonlinearity of the solved complicated equation into a treatable kind of nonlinearity. Next, we discuss the use of composite functions in the steps of the algorithms. Special attention is given to the role of the simple equation in the SEsM. The connection of the methodology with other methods for obtaining exact multisoliton solutions of nonlinear differential equations is discussed. These methods are the Inverse Scattering Transform method and the Hirota method. Numerous examples of the application of the SEsM for obtaining exact solutions of nonlinear differential equations are demonstrated. One of the examples is connected to the exact solution of an equation that occurs in the SIR model of epidemic spreading. The solution of this equation can be used for modeling epidemic waves, for example, COVID-19 epidemic waves. Other examples of the application of the SEsM methodology are connected to the use of the differential equation of Bernoulli and Riccati as simple equations for obtaining exact solutions of more complicated nonlinear differential equations. The SEsM leads to a definition of a specific special function through a simple equation containing polynomial nonlinearities. The special function contains specific cases of numerous well-known functions such as the trigonometric and hyperbolic functions and the elliptic functions of Jacobi, Weierstrass, etc. Among the examples are the solutions of the differential equations of Fisher, equation of Burgers–Huxley, generalized equation of Camassa–Holm, generalized equation of Swift–Hohenberg, generalized Rayleigh equation, etc. Finally, we discuss the connection between the SEsM and the other methods for obtaining exact solutions of nonintegrable nonlinear differential equations. We present a conjecture about the relationship of the SEsM with these methods.


Short Introduction and Overview of the Literature
Nature and human societies are full of complex systems. Atomic chains and lattices, biological systems, the dynamics of research groups, networks, social sciences, and economics are some of the many examples [1][2][3][4][5][6][7][8]. These complex systems require sophisticated research instruments and powerful computers for their study. The research on complex systems is complicated as many of these systems are nonlinear. Examples of this nonlinearity can be seen in fluid mechanics, solid-state physics, etc. [9][10][11][12][13]. Time-series analysis and models based on differential or difference equations are often used in the studies of nonlinear complex systems [14][15][16][17][18][19].
Much attention has been devoted to the research of nonlinear waves in dispersive systems [61][62][63][64][65][66]. These waves often belong to a special class of nonlinear waves: the solitons . Solitons are predicted and observed in many systems such as conducting polymers [97], plasmas [98,99], and helium films [100], and can be used, for example, in long-distance communications [101].
Below, we describe a methodology for obtaining exact solutions of nonlinear differential equations. The methodology is called the SEsM (Simple Equations Method). The methodology is a recent development in a long line of studies of methods for obtaining exact solutions of nonlinear differential equations. The classic approach is the Inverse Scattering Transform method. The basis for the Inverse Scattering Transform method was described in the paper by Gardner, Greene, Kruskal, and Miura [102]. An excellent introduction to the methodology is [103,104]. Additional results on the methodology can be found in [105]. Lax [106] proposed a general principle for associating the evolution of nonlinear equations with linear operators so that the eigenvalues of the linear operator are integrals of the nonlinear equation. We also note the contribution of Zakharov and Shabat [107,108] and the results of Fokas and Ablowitz [109,110]; Manakov [111,112]; Dod and Bullough [113]; Ablowitz, Kruskal, and Segur [114]; Kaup and Newell [115,116]; Beals [117]; Kaup [118,119]; and other authors [120,121]. This short list of contributors is in no way complete.
The Inverse Scattering Transform method stimulated much research on nonlinear differential equations. We have already mentioned above the nonlinear waves in dispersive systems and solitons. Additional results have been obtained on the integrability of dynamical systems [122,123] and the integrability of evolution equations [124][125][126][127][128][129][130][131], the construction of higher-dimensional integrable systems [132], the spectral transformations method for solving nonlinear evolution equations [133], etc. [134,135]. Other methods for obtaining exact solutions of nonlinear partial differential equations have also appeared. A prominent example of such a method is the method of Hirota [136][137][138][139][140][141].
Numerical analysis of nonlinear waves is of special interest [419]. Several results in this research area are connected to methods for the symplectic integration of Hamiltonian systems [420], the finite difference and finite element methods for the Korteweg-de Vries equation [421][422][423][424], the Galerkin method [425][426][427][428][429], and the collocation and radial basis function methods for the Korteweg-de Vries equation [430]. The numerical methods have been successfully applied to, for example, the study of the evolution of solitary waves over a shelf [431], Lax pairs symbolic computations [432], the nonisospectral scattering problems [433], the numerical solution of the nonlinear Schrödinger equation [434], the numerical study of internal wave solitons [435], and the study of the soliton decay of nonlinear Alvfen waves [436].
This concludes the brief survey of the literature, noting that it should be considered an introduction to the vast area of research on nonlinear differential equations and the methods for obtaining exact analytical solutions of such equations. Next, we present a brief history of the research that led to the appearance of the methodology, which is the focus of this review article.
As we have seen, research on the methodology for obtaining exact solutions of nonlinear differential equations is important. In the early years of these studies, there was an attempt to remove the nonlinearity of the studied equation by an appropriate transformation. An example is the Hopf-Cole transformation [437,438]. This transformation converts the nonlinear Burgers equation to the linear heat equation. Another transformation connected the nonlinear Korteweg-de Vries equation to the linear Schrödinger equation. Efforts to obtain such transformations led to the method of the Inverse Scattering Transform [439], as well as the method of Hirota [136,141]. The truncated Painlevé expansions also supplied suitable transformations of the nonlinearities [440][441][442][443]. The possibilities connected to the Painlevé expansions were also explored by Kydryashov [444], who used a truncated Painlevé expansion. In these expansions, the truncation occurs after the constant whose term is kept in the expansion. Kudryashov [445] formulated the Method of the Simplest Equation (MSE). The method is based on the determination of the singularity order n of the solved NPDE. Then, a particular solution of this equation is searched as a power series of solutions of a simpler equation called the simplest equation. For further results on this methodology and its applications, see [446][447][448][449][450][451][452][453][454][455][456][457][458][459].
About 13 years ago we started the development of a methodology for obtaining exact and approximate solutions of nonlinear differential equations. This methodology is now called the SEsM (Simple Equations Method) [460][461][462][463][464][465][466][467][468][469][470][471]. Elements of this methodology were used in our articles written decades ago [472][473][474][475]. In 2009 and 2010 [476,477], we used the ODE of Bernoulli as the simplest equation [478] in the first version of the method, the Modified Method of the Simplest Equation (MMSE), and applied the methodology to population dynamics and ecology [479]. The MMSE [480,481] uses the concept of balance equations for the fixation of the simplest equation. After this, the searched solution is constructed as a truncated power series of the solution of the simplest equation. This methodology led to equivalent results with respect to the Method of the Simplest Equation mentioned above.
The MMSE was actively used by us until 2018 [482][483][484][485][486][487]. A special role in this period was played by article [486]. There, the methodology of the MMSE was extended by the use of a class of equations d k g dξ k l = m ∑ j=0 d j g j as the simplest equations. Here, k = 1, . . . , l = 1, . . . , and m and d j are parameters. The solution of this equation contains specific cases of many well-known functions such as hyperbolic functions, trigonometric functions, the elliptic functions of Weierstrass and Jacobi, etc. The capacity of the methodology has been extended over the years. The current version of the methodology (SEsM) can use more than one simple equation. The SEsM based on two simple equations was used in [488]. The first discussion of the SEsM was in [460]. Further discussions on the SEsM can be seen in [461][462][463][464][465]. Applications of specific cases of the SEsM can be seen in [489][490][491][492][493][494][495].
Below, we discuss the SEsM in detail. The discussion covers the use of transformations to convert the nonlinearity of the solved equation to a more treatable kind of nonlinearity, the use of composite functions in the SEsM and the application of the Faa di Bruno formula for derivatives of these functions, the role of simple equations and a special class of simple equations that define an interesting special function, as well as a series of polynomials, which appear frequently in the application of the SEsM.
The text is organized as follows. In Section 2, the methodology of the SEsM is described. In Section 3, we discuss the use of the transformations of the nonlinearity of the solved equation using the SEsM. These transformations convert various kinds of nonlinearities to a polynomial kind of nonlinearity which is more easily treatable by the SEsM. In Section 4, the use of the composite functions of the SEsM is discussed with examples of such use. In Section 5, the role of simple equations in the SEsM is described. In Section 6, the connection of the SEsM to the inverse scattering transform method and the method of Hirota is discussed. Section 7 is devoted to the specific application of the SEsM for obtaining the exact solution of a differential equation connected to the SIR model of epidemic spreading, with a possible application for analysis of epidemic waves of a class of epidemics, including COVID-19. In Section 8, additional applications of the SEsM are presented. In Section 9, we discuss the conjecture about the connection of the SEsM with other methods for obtaining specific solutions on nonlinear nonintegrable differential equations. Several concluding remarks are summarized in Section 10.

The Simple Equations Method (SEsM)
The SEsM has the goal of obtaining exact and approximate solutions of nonlinear differential equations. The algorithm of the SEsM is designed for solving systems of n nonlinear differential equations. The solution is built by known exact solutions of m simpler differential equations [470]. Here, we describe a specific case of the SEsM for the solution of a single nonlinear differential equation based on the solution of m simpler equations. The MMSE is the specific case. There, m = 1.
The discussed version of the SEsM has four steps. The solved differential equation is A[u 1 (x, . . . , t), . . . , u n (x, . . . , t), . . . ] depends on the functions u 1 (x, . . . , t), . . . , u n (x, . . . , t) and some of their derivatives. The functions u i can depend on several spatial coordinates. The idea of the SEsM is to transform (1) to Here, B i are some functions. r i are the relationships among the parameters of the solved equation and the parameters of the solution. We set Thus, a system of nonlinear algebraic equations is obtained. Each nontrivial solution of (3) leads to a solution of the solved differential equation.
The steps for the transition from (1) to (3) are as follows. First, the transformation is applied. T(F, G, . . . ) is a composite function of the other functions F, G, . . . . These functions can be the functions of several spatial variables and time. If possible, the transformation can remove the nonlinearity of the solved equation. An example is the Hopf-Cole transformation, which reduces the nonlinear Burgers equation to the linear heat equation [437,438]. The removal of the nonlinearity is rarely achieved. In most cases, the transformation T transforms the nonlinearity of the solved equation to an easily treatable kind of nonlinearity such as polynomial nonlinearity. Examples of the transformation T are u(x, t) = 4 tan −1 [F(x, t)] for the sine-Gordon equation and u(x, t) = 4 tanh −1 [F(x, t)] for the Poisson-Boltzmann equation [472][473][474][475]. The Painlevé expansion can be another appropriate transformation. u(x, ..., t) = F(x, ..., t) (no transformation) is also a possibility for certain classes of nonlinear differential equations. The application of (4) to (1) leads to a treatable nonlinear differential equation for F i , G i , . . . . The transformation T may remain unfixed in the first step of the SEsM. Thus, the function T will be unknown and we have to determine it in some of the other three steps of the methodology.
Step 2 of the SEsM deals with the selection of the functions F(x, ..., t), G(x, . . . , t), . . . . They are the composite functions of the known solutions of simpler differential equations. We can let the form of the composite functions be undetermined. This, together with the form of the transformation T, can be determined in Step 3 of the SEsM. Often, F, G, . . . , are fixed in Step 2 of the SEsM. A form that is often used is . . are parameters. The relationship used by Hirota [141] is a specific case of (5). Another specific case of (5) is the power series F = N ∑ i=0 µ n f n (µ is a parameter).
This is frequently used in the MMSE.
Step 3 of the SEsM is often the most important in the application of the methodology. Here, one can determine the forms of the more simple equations and the composite functions T, F, G, . . . . All these are determined from the requirements for the satisfaction of the balance equations. The balance equations occur due to the requirement that any of r i from (3) must have at least two terms in order to ensure that a nontrivial solution of the solved equation is obtained. This determines the shapes of all composite functions and the forms of all used simpler equations.
The satisfaction of the balance equations in Step 3 leads to a system of nonlinear algebraic equations (3). This system contains relationships between the parameters of the solved equation and the parameters of the solution. In Step 4 of the SEsM, one solves the system (3). Each nontrivial solution of (3) leads to a nontrivial solution of the solved Equation (1).

The SEsM and the Transformation of the Nonlinearity of the Solved Equation
Step 1 of the SEsM is connected to the transformation of the nonlinearity of the solved equation (if necessary). If the nonlinearity of the solved equation is a polynomial one, there is no need for a transformation. Steps 2-4 of the SEsM can be applied directly in this case. In the case of the solved equation, which contains nonpolynomial nonlinearity, in Step 1 of the SEsM, one can attempt the application of a transformation that can transform the nonpolynomial nonlinearity to a polynomial nonlinearity. The following proposition can be performed, which was proved in [469].

Proposition 1.
Let us consider a differential equation for the function u(x, . . . , t), which contains two kinds of terms:

1.
Terms containing only derivatives of u; 2.
Terms containing one or several nonpolynomial nonlinearities of the function u and these nonpolynomial nonlinearities are of the same kind.
Let u = T(F) be a transformation with the following properties: 1. Property 1: The transformation T transforms any of the nonpolynomial nonlinearities to a function that contains only polynomials of F.

2.
Property 2: The transformation T transforms the derivatives of u to terms containing only the polynomials of the derivatives of F or the polynomials of the derivatives of F multiplied or divided by the polynomials of F.
Then, the transformation T transforms the studied differential equation into a differential equation containing only the polynomial nonlinearity of F.
The above proposition allows us to transform numerous nonpolynomial nonlinearities into polynomial nonlinearities. A list of some nonpolynomial nonlinearities that can be transformed into polynomial nonlinearities can be found in [469]          One illustrative example of the application of a transformation in Step 1 of the SEsM is for the equation bu 2 xx + du 2 tt = l sin 2 (u). (6) b, d, and l are parameters. In Step 1 of the SEsM, we use the transformation u = 4 tan −1 (F). The result is the conversion of (6) to an equation that contains only polynomial nonlinearities: Step 2 of the SEsM is connected with the choice of F as a composite function of simpler functions. In our case, the choice is a specific case of the composite function F( In Step 3 of the SEsM, one has to determine the form of the functions T 1 and T 2 in such a way that the balance equation leads to a balanced system of nonlinear algebraic relationships. The choice is Above, δ i and i are parameters. The balance equation has a specific solution For simplicity, specific cases of the above simple equations are considered, namely The form of the solution and the choice of the simple equations lead to the following balance system of nonlinear relationships in Step 4 of the SEsM: A nontrivial solution of this system is v, A, l, b, d, α, γ are free parameters that satisfy the condition −b(dγ 2 v 2 − l) ≥ 0. The corresponding solution of (6) is We note that in (12), one can group the free parameters. These groups are A, αv, and . Thus, the free parameters for Solution (12) are reduced to three. Let us obtain this result by obtaining (12) in another way. Above, we used a composite function of the kind F(x, t) = AT 1 (µ)T 2 (ξ); µ = αx, ξ = γt, in other words, we searched for standing waves. We can obtain (12) based on a search for the traveling wave solutions of (6), where u(x, t) = u(z) = u(x − ct). c is the velocity of the traveling wave. We introduce the parameter a 2 = l/(b + dc 2 ) and the velocity becomes Then, (6) can be factorized: u 2 zz − a 2 sin 2 (u) = (u zz + sin(u))(u zz − sin u ) = 0, and this leads to the following solution of (6): In (14), sn and dn are the Jacobi elliptic functions of modulus k, and z 0 is a constant of integration. Now, let k = 1. The Jacobi elliptic functions reduce to hyperbolic functions and (14) is reduced to its specific case Here, we have a and A = 2 exp(a 1/2 z 0 ) as free parameters. However, we have another free parameter and it is connected to the group of parameters l−ba 2 da 2 . For this group of parameters, we have the relationship (13). Let us fix a and z 0 . This does not fix c in (13) as l, b, and d are not fixed. This (15) has effectively three free parameters, similar to (12). We are grateful to the anonymous reviewer who suggested discussing the number of free parameters for the solution (12) using the considerations described in (13)-(15).

Composite Functions and their Role in the Algorithm of SEsM
One encounters composite functions in several places of the algorithm of the SEsM. The presence of composite functions leads to preferences in the choice of the simple equations and the occurrence of specific sets of polynomials in the SEsM. Below, we discuss this in detail.
The differential equations contain derivatives of unknown functions. If the unknown functions are composite, then one has to use the Faa di Bruno formula for the derivative of a composite function. The Faa di Bruno formula is written below for the function h(x 1 , . . . , x d ) of d independent variables x 1 , . . . , x d , which is a composite function of m other functions Ordering of vector indexes: for two vector indexes µ = (µ 1 , . . . , µ d ) and ν = (ν 1 , . . . , ν d ), we have µ ≺ ν when one of the following holds: In addition, the following notation is used: Based on the above notation, the Faa di Bruno for the composite derivative of a function containing the functions of many variables is [496] In (18) n =| ν | and For example, the derivatives of the composite function of two independent variables x 1 and The specific case of a composite function containing one function of one variable h = f [g(x)] is of special interest, as it is often used in practice. For this case, the Faa di Bruno formula is In (21), h (n) = d n h dx n ; f (k) = d k f dg k ; and g (i) = d i g dx i . In addition, p(n, k) = {λ 1 , λ 2 , . . . , λ n } is a set of numbers satisfying Two results of the composite functions are of interest to the methodology of the SEsM. The first is connected to the use of simple equations for exponential functions [470] Theorem 1. Let us consider a nonlinear partial differential equation that contains a polynomial P of the function h(x 1 , x 2 ) and its derivatives. The relationship for this equation is Above, N can be any natural number. We search for the solution of the above equation in the form where h is a polynomial of the functions g (1) (x 1 , x 2 ), . . . , g (m) (x 1 , x 2 ). Let each function g (i) (x 1 , x 2 ) satisfy the simple equation where α i,j is a constant parameter. Then, the solved nonlinear PDE is reduced to a polynomial of the functions g (1) (x 1 , x 2 ), . . . , g (m) (x 1 , x 2 ) containing monomials multiplied by some coefficients. These coefficients are nonlinear algebraic relations between the parameters of the solved equation and the parameters α i,j . We set the above coefficients to 0. The result is a system of nonlinear algebraic equations. Any nontrivial solution of this system leads to a solution of the solved nonlinear PDE.
The second result is connected to the use of the composite function of a function of a single variable. In this case, one can use a more complicated simple equation in the SEsM. In more detail, we consider a nonlinear partial differential equation with nonlinearities that are polynomials of the unknown function h(x, t) and its derivatives. We search for a solution of the kind where µ and ν are parameters and H is a composite function of another function g: The assumption is that f is a polynomial of g: and the general form of the simple equation is Equation (27) defines the function V a 0 ,a 1 ,...,a m (ξ; k, l, m). In this notation, k is the order of derivative of g, l is the degree of the derivative, and m is the highest degree of the polynomial of g in the defining ODE. V has as specific cases the trigonometric, hyperbolic, and elliptic functions of Weierstrass and Jacobi, etc. The result we describe is connected to a specific case of the equation for the V-function V a 0 ,a 1 ,...,a m (ξ; 1, 2, m) : The result is [486].

Theorem 2. If g 2
(1) is given by Equation (28) and f is a polynomial of g given by Equation (26), then for h[ f (g)], the following relationship holds: where K n (q, m)(g) and Z n (q, m)(g) are polynomials of the function g(ξ).
Note that for some values of n, one of the polynomials K n (q, m) or Z n (q, m) can be equal to 0.
The theorem states that the derivatives of the composite function h are constructed by specific polynomials. These polynomials are calculated by the recurrence relationships [486] with K 0 = q ∑ r=0 b r g r and Z 0 = 0. In this way, rb r g r−1 , and, for example, Another kind of simple equation leads to another set of polynomials. The simple equation (n and c j are constant parameters) is frequently used. In this case [470], one obtains the following set of polynomials. One b r q r and uses the recurrence relationship The derivative of the composite function is h (n) = L n (g), for example, Let us consider two examples of the application of the SEsM with the participation of composite functions. The first example is for the generalized Korteweg-de Vries equation In (34), p is a positive integer number and A is a parameter. The solution of Equation (34) is searched as a composite function u = h[ f (g(ξ))], where ξ = αx + βt, g is the solution of the simplest equation (28) and f is given by Equation (26). The substitution of u = h[ f (g(ξ))] in Equation (34) and the application of the SEsM leads to an equation of the kind with

and the first equation of (35) is satisfied. The second equation needs balance.
It is m = 2 + pq. For the case q = 1, m = 2 + p. Then, from Equations (26) and (28) one a j g j . Thus, the second equation from (36) is reduced to the system of nonlinear algebraic relationships among the parameters of Equation (34) and the parameters of the solution (Step 4 of the SEsM): δ is the delta-symbol of Kronecker.
The solution of Equation (37) is obtained for b 0 = 0. Then, one obtains the system of algebraic equations: The solution is . This fixes the simple equation, Equation (28), as The solution of (39) is The solution of (34) becomes In a similar manner, one can obtain numerous solutions of the Olver equation [486] ∂u ∂t

The Role of The Simple Equations in The SEsM
In general, the form of the solution of the solved nonlinear differential equation is connected to the form of the solution of the used simple equations by means of the balance equations, which have to be satisfied in Step 3 of the SEsM. We are not obligated to choose the form of the simple equations before Step 3 of the SEsM. However, if we do, additional restrictions are imposed on the form of the searched solution of the solved equation [483]. Below, we discuss this in more detail.
We consider a simple equation of the kind Using the solution Φ of this simple equation, we build a function F = F(Φ), which has to solve the more complicated equation We can choose Q and F. Thus, we have a function and the question is as follows: What is the manifold of equations that have this function as a solution? Below, are the responses to this question for the case of simple equations of the kind where and β are positive integers and the function F(Φ) is a polynomial For the construction of equations of the kind (43), one needs to know the derivatives of F with respect to ξ. Some derivatives are Derivatives (46) which has the solution In (52), a 2 Φ(ξ) 2 < d 2 and ξ 0 is a constant of integration. Equation (52) can be used for the construction of traveling wave solutions of numerous nonlinear PDEs. Solution (51) leads to specific forms of Derivatives (46)- (50). The simplest specific case of Equation (45) is One obtains Derivatives (46)-(50): The function (53) is a solution of numerous nonlinear PDEs. Let us demonstrate the construction of some of these equations. The first equation is For a solution of the kind (53), the derivatives of F are given by Equations (54) and (55). The substitution of (53)- (55) in Equation (59) leads to the proof of the following proposition.

Proposition 2 ([483]
). Each solution of the system of nonlinear algebraic relationships leads to an exact solution of the kind (53) of Equation (59).
Let us set, for example, q 3 = q 4 = q 6 = 0. Then, we obtain the conditions of a function of the kind (53) to be the solution of the equation For this, the system of the following nonlinear algebraic equations must be satisfied Let us now set q 6 = 0. Thus, we search for a solution of One solution of (60) is The substitution of (64) in (53) leads to the following solution of Equation (63) The solution (65) is an exact solution of each nonlinear partial differential equation that can be reduced to Equation (63).
One can consider more complicated differential equations such as and to search for conditions for the existence of a solution of the kind (53). One can proceed by using more complicated forms of F(Φ). One such form is dF dΦ = p 1 + 2p 2 Φ. Then, (52). The substitution of (52) and (67) in (66) leads to a system of seven nonlinear algebraic equations that have to be solved in order to construct a solution of the kind (67) and (52) of (66). A nontrivial solution of this algebraic system is Thus, the equation has the exact solution All nonlinear partial differential equations that can be reduced to (69) have the solution (70). We can continue the process by taking, for example, the equation of Bernoulli as the simple equation. The interested reader is directed to [483] for more details on this.

The SEsM and the Method of Hirota and Inverse Scattering Transform Method
The SEsM can be connected to other famous methods for obtaining exact solutions of nonlinear partial differential equations such as the method of Hirota and the Inverse Scattering Transform (IST) method. Below, we discuss the connection between the SEsM and the method of Hirota and between the IST and the SEsM for the cases of the Kortewegde Vries equation and the nonlinear Scrödinger equation.

The SEsM and the Method of Hirota
The Hirota method [136][137][138][139][140][141] is a simple and popular method for obtaining multisoliton solutions of nonlinear PDEs, which are integrable. Usually, in Step 1 of this method, a transformation of the solved NPDE is made. Frequently, the form of this transformation is In Step 2 of the method, one searches for a solution in the form where α and are parameters. Solution (72)  Hirota introduced the bilinear operators, which can be written as [497] D m These operators allow bilinear operators to have useful properties such as finding a solution for the sequence of equations for the powers of in an effective way. The Hirota method is connected to the SEsM in the following way [460]: One example of the demonstration of the connection between the SEsM and the method of Hirota is in obtaining the three-soliton solution of the Korteweg-de Vries equation: Step 1 of the SEsM is to transform the nonlinearity of this equation using (71). We obtain which can be written as (D x D t + D 4 x )( f · f ) = 0.
Step 2 of the SEsM is to fix the composite function for the solution of the equation. This is made by (72). The substitution of (72) in (75) leads to the differential equations for f 1 , f 2 , . . . : etc. The next step of the SEsM is to choose the simple equations for f 1 , f 2 , . . . . These are equations for exponential functions. One starts with dg 1 dη 1 = g 1 ; g 1 = exp(η 1 ), where η 1 = λ 1 x + ω 1 t + σ 1 and λ 1 , ω 1 and σ 1 are parameters and chooses f 1 as The substitution of (79) in (76) leads to the dispersion relation ω 1 + λ 3 1 = 0. f 1 satisfies (77) and f 2 can be taken to be 0. f 3 , f 4 , . . . can also be set to 0. One obtains The parameter is absorbed by σ 1 and the solution for f corresponds to the single soliton solution for u by the use of (71).
For obtaining the two-soliton solution of the Korteweg-de Vries equation, one uses the solution of (76) constructed by the solutions of two simple equations dg 1 dη 1 = g 1 ; g 1 = exp(η 1 ); dg 2 dη 2 = g 2 ; g 2 = exp(η 2 ). Here, η i = λ i x + ω i t + σ i and λ i , ω i and σ i , i = 1, 2 are parameters. Then, Two algebraic equations occur from the substitution of (81) in (76): The substitution of (81) in (77) leads to an equation for f 2 , for which the solution is , can be set to zero and one arrives at the two-soliton solution of the KdV equation Here, and 2 are absorbed by σ 1,2 and a 1,2 .
For obtaining the three-soliton solution of the Korteweg-de Vries equation, one needs a solution for f 1 constructed using the solutions of three simple equations for the exponential functions Here, η i = λ i x + ω i t + σ i and λ i , ω i and σ i , i = 1, 2, 3 are parameters. f 1 is the sum of the solutions of the three simple equations The substitution of (83) in (76) leads to three dispersion relations ω i + λ 3 i = 0, i = 1, 2, 3. For f 2 , after a substitution of (83) in (77), one obtains f 2 = a 12 exp(η 1 + η 2 ) + a 13 exp(η 1 + η 3 ) The substitution of f 1 and f 2 in (78) leads to the following relationship for f 3 : 23 . This leads to f 4 = 0, f 5 = f 6 = · · · = 0. The solution for f is and this leads to the three-soliton solution of the Korteweg-de Vries equation.
We skip Step 1 of the SEsM (the transformation of the nonlinearity). In Step 2 of the SEsM, we present u(x, t) as a composite function of other functions u 1 , u 2 , . . .
Initially, one treats as a small parameter. The substitution of (86) in (85) leads to equations for u 1 , u 2 , . . . as follows: ∂u n ∂t and the equation for u 1 is ∂u 1 ∂t The solutions of the obtained equations for u 1 , u 2 , . . . are connected through a Fourier series to the solution of the simple equation for exponential functions. For u 1 we use where dλ(k) is an appropriate measure in the complex plane C and the term (−k) is introduced for convenience. The form of u 1 influences the relationships for the equations for u 2 , u 3 , . . . . One obtains The solutions for u 2 , u 3 , . . . are assumed to be of a form similar to that of u 1 : The substitution of (92) in the equations for u 2 , u 3 , . . . leads to the fixation of Φ n and Ω n . The fixation leads to ; ; Ω 4 = (k 1 + k 2 + k 3 + k 4 )x + (k 2 1 + k 2 2 + k 2 3 + k 2 4 )t. (93) This leads to the following relationship for u n and the solution of the Korteweg-de Vries equation is One rewrites (95) as ; ω(k) = −k 3 . Two possibilities exist for the measure dλ(k), discrete or continuous. For the case of discrete dλ(k), the integral can be replaced by a sum: Here, P mq = p m p q k m +k q ; p m = a mp (ik m ) = a m exp 1 2 (−k m x + k 3 m t) , and (97) can be written as where p is the column vector of all p m , P is the square matrix of all P mq , and p T is the transpose of p. The matrix P is real symmetric and positive definite for k m > 0 and real for a m . Thus, u is nonsingular for > 0 and can be absorbed into the coefficients a m . Because of this, Solution (98) is not limited to only small values of . Solution (98) can be rewritten in a more common form. Observing that ∂P ∂x = − 1 2 pp T , we arrive at This is the relationship for the multisoliton solution of the Korteweg-de Vries equation.
For the case of continuous measure dλ(k), we obtain the Gelfand-Levitan-Marchenko equation from the methodology of the Inverse Scattering Transform method as follows. P is considered more a general operator than a square matrix: (P f )(k) = C dλ(l)P(k, l) f (l).

Thus, the solution of the KdV equation is
Solution (100) can be written as Here, K(x, y) = − 2 p T (x)[I + P(x)] −1 p(y), and after a calculation, one obtains This is the GLM equation with

The SEsM and the Inverse Scattering Transform Method: The Case of the Nonlinear Schrödinger Equation
We consider the nonlinear Schrödinger equation Step 1 of the SEsM (the transformation of the nonlinearity) is skipped. In Step 2 of the SEsM, we present ϕ(x, t) as a composite function similar to (86). Performing the calculations for the solutions of the equations for the components of ϕ using differential equations for exponential functions as simple equations (Steps 3 and 4 of the SEsM) we arrive at In (105), Ω n n ∑ 1 [k j x + (−1) j k 2 j t]. dλ(k) and dµ(k) are measures such that dν * (−k * ) = dλ(k) ( * denotes complex conjugation). dλ j = dλ(k j ); dµ j = dµ(k j ). There is a quantity analogous to P from the case of the Korteweg-de Vries equation above. The definition of this quantity for the case of the nonlinear Schrödinger equation iŝ In Ω n , one has two dispersion relationships: ω 1 (k) = k 2 and ω 2 (k) = −k 2 . Inp, we take the first of them. Denoting the conjugate operator of P as P, one can rewrite (105) as For | ϕ | 2 , one obtains This is reduced to the relationship of Zakharov and Shabat [500] for the multiple-envelope solutions of the nonlinear Schrödinger equation for the specific case of (108) of discrete measure (discrete λ(k)) and when p is a vector and P is a square matrix. In this case,

Special Application: The SEsM and the SIR Model of Epidemics
Below we discuss the exact solutions of an equation. This equation is obtained based on the SIR model in epidemiology. The obtained solutions can describe epidemic waves caused by different diseases (COVID-19-inclusive).
The equation is obtained using the classic idea of Kermak and McKendrick [501] for the transformation of the SIR model with constant coefficients to a single nonlinear differential equation. We consider an epidemic in a population. The population is divided into three groups: susceptible individuals-S; infected individuals-I; recovered individuals-R. The model equations for the time changes in the numbers of individuals in the above three groups are In (refsir1), τ is the transmission rate and ρ is the recovery rate. These rates are assumed to be constants.
is the total population, which is assumed to be constant. The model (110) can be reduced to a single equation for R as follows. From the last equation of (110) we have The substitution of (112) in the first equation of (110) leads to the relationship Here, S(0) and R(0) are the numbers of susceptible and recovered individuals at the time t = 0. The substitution of (111) and (113) in the last equation of (110) leads to the differential equation Below, we assume R(0) = 0 (no recovered individuals at t = 0). Let us consider the ratio τR ρN . The maximum of the ratio R/N is 1 and if τ << ρ, then τR ρN << 1. Let us consider a kind of epidemic, where τ << ρ (transmission rate is much lower than the recovery rate). Then, exp − τ ρN R can be represented as a Taylor series M has infinite value in the full Taylor series but we can truncate it at We apply the most simple version of the SEsM to (116). The solution is searched as a composite function of a single function U and the form of the composite function is a polynomial where µ l and K are parameters. We also fix the form of the simple equation for U where ν l and L are parameters. In order to proceed with the application of the SEsM, we have to obtain a balance equation that connects the parameters M, K, and L. The form of this balance equation is Then, we proceed as follows. We fix the value of M. This fixes the form of (116). Then, we fix the value of K. This fixes the form of the composite function for the solution of the equation. Then, from (119) we obtain the value of L. The fixation of the parameters M, K, and L and the substitution of (117) and (118) in (116) lead us to the system of nonlinear algebraic equations as required by the algorithm of the SEsM. The nontrivial solutions of this system lead to solutions of R. The substitution of these solutions in (112) leads to the corresponding solutions for I. Finally, the solutions for R and I are substituted in (111) and we obtain the solution for S.
Below, we consider the simplest case M = 2. This is the classical Kermak-McKendrick case. From (116), the equation we have to solve is where (120) is an equation of the Riccati kind. From (119), Let us consider the case K = 1, L = 2. In this case, and the simple equation is The substitution of (122) and (123) in (120) leads to a system of nonlinear algebraic equations. One nontrivial solution of this system is µ 0 = 0, µ 1 = 1; ν 0 = γ; ν 1 = β, ν 2 = α. The simple equation, Equation (123), coincides with (120). As we consider the case τ << ρ then, A particular solution of (120) is where θ 2 = β 2 − 4αγ > 0 and C is a constant of integration. The substitution of (124) in (112) leads to the well-known bell-shaped curve for the evolution of the number of infected individuals I Note that − θ 2 2αρ is positive as α is negative. We know a particular solution (124) of (120). Then, we can write the general solution of (120) as where D is a constant and v(t) is the solution of the linear differential equation The solution of (126) is where E is a constant of integration. Then, the general solution of Equation (120) is This solution introduces corrections to the well-known curve for the time evolution of recovered persons within the scope of the SIR model. The substitution of (128) in (112) leads to the following relationship for the time evolution of the number of infected persons: This solution leads to corrections to the well-known bell-shaped curve for the time evolution of infected individuals in the SIR model. For example, let t = −C. Then, from (125), and from (129), The presence of the general solution (128) of (116) makes the further consideration of the subcases of the case M = 2 unnecessary. What we can obtain by considering these cases are the specific cases of the general solution (128).

Other Examples of the Application of the SEsM
The SEsM and its most simple version, the MMSE, have numerous applications. Below, we mention several results. In most cases, we search for traveling wave solutions of the kind u(x, t) = u(ξ) = u(x − vt), and the composite function from Step 2 of the SEsM is given by the power series of the solution Vof the simple equation In (132), ν > 0, µ > 0, θ µ are parameters. In most cases, we use as the simple equation where γ α is a parameter. Equation (133) where k is an integer and k > 1. The solutions that are used are The solution that is used is Here, θ 2 = b 2 − 4ac > 0 and ξ 0 is a constant of integration. As a specific case of the use of the equation of Riccati, we also consider the extended tanh-function equation Equation (138) is obtained from (136) when b = 0, a = −a 2 , and c = c 2 . The solution of (138) that is used is where a 2 V(ξ) 2 < c 2 and ξ 0 is a constant of integration. The first illustration of the methodology of the SEsM is for obtaining traveling wave solutions of the class of equations In (140), α p , β q , and µ m are parameters. Specific cases of (140) are used as model equations in the theory of the migration of populations.
Such kinds of equations are widely used in the modeling of genetic diffusion.
One searches for traveling wave solutions Q(x, t) = Q(ξ) = Q(x − vt). Taking into account that N = max(N 1 , N 2 ), one can extend (140) to where the coefficients of the additional terms are equal to 0. The traveling wave solution and ν n = β n + α n (−v) n reduce (141) to As the simple equation, we use the equation of Bernoulli (134). The substitution of (132) and (134) in (142) leads to a polynomial P = κ r φ r + κ r−1 φ r−1 + · · · + κ 0 = 0. In P, r is an integer. The coefficients κ depend on the parameters of the solved equation and the parameters of the solution. The solution is obtained when all coefficients are set to 0. The result is a system of nonlinear algebraic relationships Next, we move to Step 3 of the SEsM, ensuring that each equation from (143) has at least two terms by balancing the highest powers arising in (142). This leads to the following balance equation (for the case when the Bernoulli equation is used as the simple equation): LM = L + N(k − 1). Thus, the choice of the highest power L of the solution polynomial leads to the fixation of the parameter k for the Bernoulli equation If the Riccati equation is used as the simple equation, then the balance equation is Several examples of exact solutions follow for cases where the Bernoulli equation is used as the simple equation. We first set N = 2, M = 2, and (142) becomes A specific case of this equation is the equation studied by Ablowitz and Zeppetela [502] Another solution of the last equation is obtained using the tanh-method. This solution is [503] We search for a solution of (145) of the kind Q(ξ) = ∑ L i=0 a i φ i ; dφ dξ = bφ (1+L/2) + aφ. L = 2 is assumed and for this case, in Step 4 of the SEsM one obtains the system of nonlinear algebraic relationships between the parameters of the solutions and the parameters of (142) A solution of (146) is a 0 = 6 25 . The corresponding solution of (145) is for the case when b < 0 and a > 0. Another solution of (146) is a 0 = a 2 1 4a 2 , ν 2 = 5µ 1 6a , . The corresponding solution of (142) is Next, we set L = 4 ( L = 3 is impossible as k has to be an integer). In Step 3 of the SEsM, we obtain the system of nonlinear algebraic relationships between the parameters 24ν 1 a 4 b 2 + µ 2 a 2 4 = 0, 2µ 2 a 4 a 3 + 15ν 1 a 3 b 2 = 0, 4ν 2 a 4 b + 40ν 1 a 4 ba + 8ν 1 a 2 b 2 + 2µ 2 a 4 a 2 + µ 2 a 2 3 = 0, 24ν 1 a 3 ba + 2µ 2 a 4 a 1 + 3ν 1 a 1 b 2 + 3ν 2 a 3 b + 2µ 2 a 3 a 2 = 0, 2µ 2 a 4 a 0 + 16ν 1 a 4 a 2 + 4ν 2 a 4 a + µ 1 a 4 + 2µ 2 a 3 a 1 + µ 2 a 2 2 + 2ν 2 a 2 b + 12ν 1 a 2 ba = 0, µ 1 a 3 + 9ν 1 a 3 a 2 + 3ν 2 a 3 a + ν 2 a 1 b + 2µ 2 a 2 a 1 + 2µ 2 a 3 a 0 + 4ν 1 a 1 ba = 0, 4ν 1 a 2 a 2 + µ 1 a 2 + 2µ 2 a 2 a 0 + 2ν 2 a 2 a + µ 2 a 2 1 = 0, ν 2 a 1 a + ν 1 a 1 a 2 + 2µ 2 a 1 a 0 + µ 1 a 1 = 0, One solution of this system is a 0 = 6 . The corresponding solution of (145) (b > 0, a < 0) is Another solution of (149) is a 4 = ba 2 2a , a 3 = 0, a 1 = 0, a 0 = aa 2 2b , µ 1 = 12 5 aν 2 , µ 2 = − 24 5 bν 2 a 2 , ν 1 = ν 2 10a . The corresponding solution of (145) is Here, b > 0, a < 0. Next, one can set N = 2, M = 3. The equation that is to be solved is This equation is connected to the Huxley equation for gene propagation [504]. Other equations having polynomial nonlinearity of the third order are the Burgers-Huxley equation and equations that model aspects of neuronal activity [505]. Several solutions of (152) are below. For the case L = 1, the SEsM leads to a solution of the kind For the case b > 0, a < 0, the solution is For L = 2, the solution of (152) for the case a > 0, b < 0 is One can continue by assuming N = 2, M = 4. The equation that has to be solved is Equations of this kind arise in the theory of the migration of populations [476,477]. The application of the SEsM leads to the conclusion that this equation has solutions of the form for the case b > 0 and a < 0. Next, we can set N = 2, M = 5. The equation that has to be solved is The SEsM leads to the conclusion that this equation has solutions of the kind Q(ξ) = ∑ L i=0 a i [φ(ξ)] i , dφ dξ = bφ (1+2L) + aφ. For the case L = 1 and for a > 0, b < 0, the solution is [478] where The solutions above are for the use of the equation of Bernoulli as the simple equation in the SEsM. One can also use other equations as simple equations. One example is the equation of Riccati. Numerous solutions of the discussed equations obtained based on the application of the SEsM based on the Riccati equation as the simple equation can be seen in [478].
The SEsM has been applied to obtain exact solutions of more complicated classes of equations. One such class is connected to modeling PDEs from ecology and population dynamics [480] A p (Q), B r (Q), C s (Q), D u (Q), and F(Q) are polynomials of Q. Specific cases of (160) have been obtained as model equations in the theory of the migration of populations [476,477]. We search for traveling wave solutions Q(x, t) = Q(ξ) = Q(x − vt). This transforms (160) into In (161), . A p , B r , and F are assumed as follows: Here, q i , q * i , and q i are coefficients. We discuss obtaining the balance equations for several specific cases below. The corresponding exact solutions can be seen in [480].
We discuss the case of the application of the SEsM for obtaining the exact solution of (161) when the simple equation is the equation of Bernoulli and the solution is searched in the form φ(ξ) is a solution of the equation of Bernoulli, which is used here as the simple equation. Thus, we fix the form of the composite function for the solution and we fix the simple equation. Next, we have to go to Step 3 of the SEsM to obtain the balance equations. In order to obtain the balance, we have to estimate the maximum powers of the term of the solved equation for the case when the equation of Bernoulli is used as the simple equation. We name the maximum powers BL 1p , BL 2r , and BL 3 . The estimations are as follows: LM 1p is the maximum power of φ(ξ) arising from the polynomial A p (Q).
We denote as max(BL 1p ) the largest value among all BL 1p and with max(BL 2r ), the largest value among all BL 2r . The theoretically possible balance equations are Often, one of these three terms is smaller than the other two. Several examples of the determination of the balance equations are below. For the equation where D is the constant coefficient of diffusion, Ψ(Q) is a polynomial of Q as follows The largest powers of φ(ξ) are BL 11 = L + k − 1; BL 12 = LM 1 + 2(k − 1); BL 22 = LM 1 + 2(k − 1).
Equation (164) has to be a nonlinear PDE. Then, M 1 > 2. This leads to BL 12 > BL 11 and BL 22 > BL 11 . The balance equation becomes The next example is the balance equation for the reaction-diffusion equation with a diffusion coefficient that depends on population density where we assume D(Q) = ∑ For the case of the use of the Bernoulli equation as the simple equation, the largest powers of φ(ξ) arising from the four terms of (166) are BL 11 = L + k − 1; BL 12 = LM 1 + L + 2(k − 1), BL 22 = L(M 1 − 1) + 2(L + k − 1); BL 3 = LM. BL 11 is smaller than BL 12 = BL 22 . Two possibilities remain. If M is such that BL 3 is smaller than BL 12 , then the balance equation is BL 12 = BL 22 , and no requirements are imposed on L and M 1 . The second possibility is to balance BL 12 and BL 3 . Thus, BL 12 The balance equations for the other differential equations of the studied class, as well as the corresponding exact solutions, can be seen in [480].
A more complicated example is connected to the extended Camassa-Holm equation [484]. The Camassa-Holm equation is a model equation for the x− component u of the fluid velocity at a certain depth z 0 below the fluid surface. The Camassa-Holm equation was introduced as a completely integrable bi-Hamiltonian dispersive shallow water equation [506]. The Camassa-Holm equation is [484,507] and δ are the amplitude parameter and the shallowness parameter and U is a quantity connected to the height of the wave η over the water surface.
Below, we present an exact solution of the generalized Camassa-Holm equation [484] ∂U ∂T p 1 ; p 2 ; p 3 ; p 4 ; p 5 are parameters. Specific cases of (169) are as below. For f (U) = 1 2 U 2 ; g(U) = C + 3U 2 (C is a constant) and p 3 = p 4 p 5 , one obtains Equation (169) contains the Camassa-Holm equation as a particular case when p 2 = 2κ; p 4 = ; p 2 = δ 2 ; p 3 = 2 δ 2 . In general, in Equation (169), f (U) and g(U) can be arbitrary polynomials of U. We skip Step 1 of the SEsM (no transformation of the nonlinearity). We search for traveling wave solutions and introduce the coordinate ξ = X − vT, where v is the velocity of the traveling wave. In Step 2 of the SEsM, we fix the composite function for the solution to be a power series of the solution of the simple equation. In Step 3 of the SEsM, we have to determine the balance equations for the corresponding simple equations. For the case when the simple equation is the equation of Bernoulli, the balance equation is where k > 1 and J > I > 1. For I > 1, another balance equation is possible. It is when (J − I)ν 1 < 2(k − 1). This balance is possible because of the existence of two terms in (169), which lead to the maximum power Iν 1 + 3(k − 1)of the terms of the resulting polynomial, and these powers are maximum powers for the entire polynomial when (J − I)ν 1 < 2(k − 1). A third balance equation is possible in the case of the use of the Bernoulli equation as the simple equation. This balance is when I = 1. The corresponding balance equation is Another possibility for the simple equation is the equation of Riccati (136). The solution of Equation (169) is, again, of the kind (132). The possible balance equations in this case are The possible balance equations for the case when the extended tanh-equation is used as the simplest equation are the same as for the case of the Riccati equation.
We return to the case of the use of the SEsM with the equation of Bernoulli as the simple equation. The values of parameters in the balance equation, Equation (171), are set to J = 3, I = 2, k = 2. Then, from (171), ν 1 = 2. The equation we solve is The solution of Equation (170) is searched in the form V(ξ) is the solution of the equation of Bernoulli for k = 2. Substituting (178) in (170), we obtain a system of seven nonlinear relationships between the parameters of the solved equation and the parameters of the solution (178).
This leads to the following solution of Equation (177): for c < 0 and a > 0. For c > 0 and a < 0, the solution is Another possibility is to use the equation of Riccati (136) as the simple equation. We use (174) for the balance equation and set the parameters to J = 1; ν 1 = 2; I = 2. The equation that is to be solved is In Step 4 of the SEsM, we obtain a nonlinear algebraic system of eight relationships between the parameters of the solved equation and the parameters of the solution. A nontrivial solution of this system is , The corresponding solution of Equation (182) based on the solution (137) of the Riccati equation is In (184), θ 2 = c 2 − 4ad > 0. Finally, we can use the extended tanh-equation, Equation (138), as the simple equation. From the corresponding balance equation, we choose the parameters I = J = 2 and ν 1 = 2. The equation to be solved is In Step 4 of the SEsM, one obtains a system of nonlinear algebraic equations that contains eight relationships. A nontrivial solution of this system is The corresponding solution of the solved equation, Equation (185), is The next example we consider is connected to the generalized Swift-Hohenberg equation. The 1+1-dimensional version is discussed: Here, r is a parameter and N(u) is a polynomial nonlinearity N(u) → ∑ s * s=0 c * s u s , where c * s , s = 0, 1, 2, . . . are parameters. Traveling wave solutions are searched. We set c 0 = c * 0 ; . . . , c s = c * s , and the solved equation becomes The balance equation that is obtained in Step 3 of the SEsM is We consider the specific case s * = 5. The corresponding Swift-Hohenberg equation becomes The balance equation is For β = 2, η = 1. β = 2 means that we use the equation of Riccati (136) or the extended tanh-equation (138) A nontrivial solution of this system is The corresponding solution of (191) is where a = −a 2 and c = c 2 .
For the case β = 3, one obtains from (192) which for γ 0 = γ 2 = 0 and γ 1 = a, as well as γ 3 = b, is reduced to the equation of Bernoulli (134). A nontrivial solution of the corresponding system of nonlinear algebraic relationships obtained in Step 4 of the SEsM is The solution of the Swift-Hohenberg Equation (191) is for the case b < 0, a > 0 and for the case b > 0, a < 0. θ 0 and θ 2 are free parameters. The next example is connected to the generalized Rayleigh equation F and G are the functions of the corresponding terms.
Equation (201) occurs in the theory of sound. We consider the case where F and G are polynomials Traveling wave solutions u(x, t) = u(ξ) = u(x − vt) are searched for the solved equation The balance equation that occurs in Step 3 of the SEsM is For the case m = 2, we solve the equation and the balance equation is Here, we set β = 2. This fixes the simple equation to be the equation of Riccati or the extended tanh-equation. In addition, we set s * = 3 and the solved equation becomes For the case of the extended tanh-equation as the simple equation, one obtains a system of nonlinear relationships in Step 4 of the SEsM. A nontrivial solution of this system is The solution of (207) is Next, we set β = 3. The simple equation, Equation (133), is again reduced to the equation of Bernoulli (134) with k = 3. Moreover, we set s * = 3. Then, from the balance equation, η = 4. The solved equation is (207), and in Step 4 of the SEsM, one arrives at the system of nonlinear algebraic relationships. A nontrivial solution of this system is The corresponding solution of the solved equation is for the case b < 0, a > 0 and for the case b > 0, a < 0. The last example illustrates the application of the SEsM to the case of the use of two simple equations. The set of equations that is discussed contains equations of the nonlinear Schrödinger kind [488]: In order to obtain a solution of Equation (213) one needs the two simplest equations. We skip Step 1 of the SEsM (no transformation of the nonlinearity, as the nonlinearity is a polynomial one). In Step 2 of the SEsM, we can choose a composite function for the solution.
q(x, t) = [V a 3 (ξ; 1, 2, 2D Here, a 3 = (a 0 , . . . , a 2D+2 ) and a 2i = c i−1 , a 2i+1 = 0, i = 0, D + 1. The above solutions can be expressed by the specific case of the V-functions for the selected values of the parameters. One example is for the solutions given by Equation (222). These solutions are:
Case B = 1. Equation (213) is The solution of Equation (225) is Here, ℘ is the elliptic function of Weierstrass. Case B = 2. Equation (213) is The solution is expressed by the Jacobi elliptic functions [488].

The SEsM and Other Methods for Obtaining Specific Exact Solutions of Nonlinear Nonintegrable Differential Equations
An interesting direction for future research is the relationship between the SEsM and other methods for obtaining exact solutions of nonlinear differential equations. We have shown above that by using the appropriate composite function of Step 2 of the SEsM and simple equations for exponential functions, the SEsM can be connected to the method of Hirota and the Inverse Scattering Transform method. These two methods are often used for obtaining exact solutions of integrable nonlinear differential equations. The situation with respect to the methods for obtaining particular solutions of nonintegrable nonlinear differential equations is as follows. It was shown that the SEsM contains as specific cases [511] the Jacobi Elliptic Function Expansion method [512], F-Expansion method [513][514][515], Modified Simple Equation method [516], Trial Function method [517,518], General Projective Riccati Equations method [519], and First Integral method [520]. Other methods that demonstrate specific cases of the SEsM are [495] the Homogeneous Balance method [521][522][523][524] and the Auxiliary Equation method [525]. The list continues [466] with the G'/G method [526,527], Exp-function method [528], tanh method [529], and the method of the Fourier series. We note that the SEsM can also deal with differential equations containing fractional derivatives [465]. The list can be continued and we present here the following conjecture.

Conjecture 1.
Any method for obtaining exact analytical solutions of nonintegrable nonlinear differential equations is a specific case of the SEsM.
The motivation for this conjecture is as follows. The SEsM is based on known solutions of simple equations and the solution of the more complicated solved equation is constructed as a composite function of these solutions. Any method for obtaining exact analytical solutions of nonlinear differential equations has to construct these solutions using known functions. These known functions are the solutions of the simple equations. The solution of the solved equation is a function of these known functions. This is the composite function used in the SEsM. Then, if one wants to find a method for obtaining the exact solution of a nonlinear differential equation that is not a specific case of the SEsM, one has to do one of the following things:

1.
Use the forms of the known functions that are not solutions of any simple equations. This will be extremely difficult as one has to use functions that are not solutions of differential equations.

2.
Construct the searched solution as a function of the known functions, which is not a composite function. This is also an extremely difficult task.

Concluding Remarks
In this review article, we discuss an effective method for obtaining exact analytical solutions of nonlinear differential equations, the Simple Equations Method (SEsM). Through a relatively large overview of the literature, we show that the research area connected to the search for exact solutions of nonlinear differential equations is of much interest because these equations are used in the mathematical models of numerous complex systems, and the mathematical aspect of the search of exact solutions is also quite interesting. The SEsM has a simple algorithm based on the appropriate transformation of the nonlinearity of the solved equation (if necessary), and then the solution is constructed as a composite function of solutions of simpler differential equations. The SEsM has four steps. In Step 1, one may use a transformation to convert the nonlinearity of the solved equation to a more treatable kind of nonlinearity. A treatable kind of nonlinearity is polynomial nonlinearity.
Step 2 of the SEsM is connected with the use of a composite function in order to build the solution of the solved equation. We show that the composite function favors the occurrence of a certain class of simple equations that define a special function, the V-function. This function contains as specific cases numerous elementary and special functions and leads to the occurrence of systems of polynomials in the process of the application of the SEsM to the solved equation. The SEsM transforms the solved equation into a sum of functions multiplied by some coefficients, which are relationships between the parameters of the solved equation and the parameters of the solution. The setting of these coefficients to 0 leads to a system of nonlinear algebraic equations. In order to obtain a nontrivial solution of the solved equation, this system of nonlinear algebraic relationships must be balanced. This leads to additional relationships called balance equations. The balance equations are important as they may lead to a fixation of the form of the simple equation and the form of the composite function for the solution of the solved equation. Thus, the first three steps of the SEsM convert the solved equation to a system of nonlinear algebraic relationships. The solution of this system is the task of Step 4 of the SEsM. Each obtained nontrivial solution of the system of nonlinear algebraic relationships leads to a nontrivial solution of the solved equation.
We note that in general, the SEsM is constructed to be able to solve systems of nonlinear differential equations. Above, we presented the version of the method for the solution of one nonlinear differential equation. We demonstrated that this version of the SEsM is connected to well-known methods for obtaining multisoliton exact analytical solutions of integrable differential equations, the Inverse Scattering Transform method, and the method of Hirota. We The reason is that we showed the application of the SEsM for the case of two simple equations in the process of obtaining exact analytical solutions of these equations.
As we have shown above, the methodology of the SEsM is effective. It is connected to well-known methods for obtaining exact multisoliton solutions of nonlinear differential equations and contains as specific cases many methods for obtaining specific exact analytical solutions of nonintegrable equations. Future research on the SEsM will be in two directions: (i) an extension of the area of the application of the methodology by treating additional classes of nonlinearities and equations; and (ii) the application of the methodology to systematic methods for obtaining exact solutions of nonlinear differential equations.
Funding: This research was partially supported by the project BG05 M2OP001-1.001-0008 "National Center for Mechatronics and Clean Technologies", funded by the Operating Program "Science and Education for Intelligent Growth" of Republic of Bulgaria.