The Functional Expansion Approach for Solving NPDEs as a Generalization of the Kudryashov and G (cid:48) / G Methods

: This paper presents the functional expansion approach as a generalized method for ﬁnding traveling wave solutions of various nonlinear partial differential equations. The approach can be seen as a combination of the Kudryashov and G (cid:48) / G solving methods. It allowed the extension of the ﬁrst method to the use of second order auxiliary equations, and, at the same time, it allowed non-standard G (cid:48) / G -solutions to be generated. The functional expansion is illustrated here on the Dodd–Bullough–Mikhailov model, using a linear second order ordinary differential equation as an auxiliary equation.


Introduction
Solving nonlinear partial differential equations (NPDEs) is an important issue in many problems from mathematical physics. This is mainly related to the fact that the integrability of these equations is a problem in itself and there are not any clear prescriptions or algorithms that can be used for solving such equations. Many approaches have been proposed, both for establishing if the equations are integrable or not, and for solving the integrable ones. It is important to mention that the same NPDE could present many classes of solutions, depending on the values of the parameters appearing in the equation.
An important class of solutions is represented by the traveling wave solutions. They are very important in the theory of solitons and are related to a symmetry transformation, which leads to a one-dimensional, nonlinear ordinary differential equation (NODE) [1]. This is obtained by using the wave variable and a whole symmetry group accepted by the initial equation. Let us consider that the variable u(t, x) defined in a 2D space satisfies an NPDE of the following form: F(u, u t , u x , u xx , u tt , · · · ) = 0. (1) Traditionally, the wave variable includes the wave velocity V, and has the following form: It transforms (1) into a NODE of the form: ∆(u, u , u , · · · ) = 0, where u = du(ξ)/dξ. In principle, solving (3) is simpler than solving (1) and then, by pulling back the solutions of (3) to the initial variables {x, t}, one can find solutions of the NPDE (1). Many approaches for finding traveling wave solutions have been proposed and are currently used in literature. Some of them have a strong theoretical basis and are related to approaches such as the inverse scattering method [2], Lax operators [3,4], Hirota and super-Hirota biliniarization [5][6][7][8], Lie symmetry theory [9,10], the ghost field method [11][12][13][14], the homotopy technique [15], etc. There are also direct approaches, trying to see if the investigated NPDE accepts traveling wave solutions with a pre-defined mathematical form: harmonic solutions expressed through sine-cosine [16], hyperbolic solutions expressed by cosh or tanh [17,18], the first integral method [19,20], etc. These attempts were generated by the fact that such solutions correspond to important equations from soliton theory, such as Riccati [21] or Jacobi [22,23] equations. These were only a step before the invention of the so-called auxiliary equation method for solving NPDEs [24]. In this case the NPDE solutions u(ξ) have to be expressed as combinations or expansions of any of the known solutions G(ξ) of these basic "auxiliary" equations. Many investigation methods based on auxiliary equations have been proposed, such as the exponential method [25,26] or the Kudryashov method [27,28]. In the last mentioned case, the following is supposed: where G(ξ) represents a solution of the Riccati equation: As the Riccati equation is a very simple, first order equation, the solution (4) depends on G(ξ) only, having at the end the following form: Many other types of auxiliary equations have been considered in the literature, some of them being of higher differential orders. For example, if a second order auxiliary equation, ∆(G, G , G ) = 0 is considered, the solution u(ξ) will also depend on the first derivative of G(ξ) and it should be expressed as follows: The question is how the Kudryashov method can be extended in the case of second order auxiliary equations. Many authors use the G /G-method, in its classic form [29], or in various extended or generalized versions [30,31]. In this approach the solution (6) has to mainly be considered as an expansion with constant or function coefficients a i of the following form: The method does not offer a clear answer as to why the only possible combination is in the form of (7) and if other extended forms are still possible: answering these questions was the main aim of this work. The starting point is represented by our previous paper [32], where the functional expansion method was proposed. This approach can be seen as an extension of many other approaches, including the Kudryashov method, and it has been shown that for models such as KdV, Gardner and Kundu-Mukherjee-Naskar (KMN) [33], more general solutions other than G /G can be generated that use a linear second order ODE as an auxiliary equation. Here, the Dodd-Bullough-Mikhailov (DBM) equation was used as an exemplifying model.
As many other methods can obtain traveling wave solutions, the functional expansion method has two important ingredients: (i) transformation of the NPDE into a NODE using the wave variable; (ii) finding solutions of the NODE in terms of the known solutions of the auxiliary equation. Both ingredients bring specific aspects and could generate intensive analysis and discussions. The NPDE solutions strongly depend on how these ingredients are chosen. We will see that the functional expansion method supposes a very specific (double) balancing procedure that creates the differences between this and the other approaches.
The paper is structured in the following sections: after these general ideas, the functional expansion method will be briefly reviewed in the next section. A second order auxiliary equation and a specific class of NPDEs will be considered, which is quite a general class of equations, including models such as Korteweg de Vries, nonlinear Schrödinger, Klein-Gordon, etc. As we already mentioned, to illustrate how the functional expansion works, a specific case belonging to the mentioned class of equations, namely the DBM equation, will be considered in the third section of the paper. We obtain traveling wave solutions for the model and we try to prove if they are new ones or if they can be reduced to already-known solutions for the same model. This checking is very important [34], and we will comment on it at the end of the paper.

The Functional Expansion Method
Let us consider that to solve (1) we transform it into a NODE of the form in (3), using the wave variable (2). We are looking for the solutions of (3) that can be expressed as combinations of solutions G(ξ) of an auxiliary equation of the form: Depending on the differentiability order of (8), the most general form for the solution of (3), should be: In particular, the previous relation could be considered as: Here P i (G) are 2m + 1 functionals depending on G(ξ) and that have to be determined. H(G, G , G , · · · ) can be a very general expression containing G(ξ) and its derivatives. Depending on the form of P and H, one can generate very complex solutions. The choices covering almost all the approaches currently used in the literature are H depending on G and G only. More strictly, H is usually considered as a formal series expansion at most linear in the two variables: The generalized Kudryashov method corresponds to the case h 0 = 0, h 1 = 0, h 2 = 0. There are other choices, too. For example, if we consider the opposite situation, h 0 = h 1 = 0 and h 2 = 0, we obtain from (10) an expression of the following form: The generalized and improved G /G method [30,31] corresponds to (12) with the following choices: The approach from [35] corresponds to (12) with P i = a i G −i + b i+1 G −i+1 , while the (w/g) method is recovered for P = 1/G and H = w(G), with an adequate choice for the auxiliary equation. The representation used in [36] is also included in (9), but it does not accept the condensed form (12). It imposes an H(G, G ) of the following form: The functional expansion method deals with solutions in the form of (12), where P i (G) are arbitrary functionals that have to be determined. This is performed by applying a balancing procedure for determining m, followed by vanishing the coefficients for different powers in G . At the end, we obtain a system of NODEs in the functionals P i (G) and its derivatives,Ṗ i ≡ dP i dG ,P i ≡ d 2 P i dG 2 , etc. If this system can be directly solved, we obtain very general solutions of the form in (12). However, it seems that direct solving is not always straightforward, so we need ways out to particular solutions. Such a way, as was proposed in [32], consists of looking for solutions given as expansions in G that can be represented as rationales with a polynomial numerator, N (G), and denominator, D(G): This choice, with the functionals {P i , i = 0, 1, · · · , m} as ratios of polynomials, gives an answer to the previous formulated question on how to generalize the Kudryashov method for second order auxiliary equations, with (15) being quite similar to (4), Kudryashov's choice. In (15), n(N i ) and n(D i ) are the degrees of N i (G) and D i (G), respectively. The parameters {π iα , ω iβ } are constants that have to be determined in order to be able to write down the effective solutions u(ξ). To determine these degrees, we need, as we will see, a supplementary balancing procedure, taking into account the "degree" attached to P i (G): In principle, these degrees can be positive, negative or zero but, considering this supplementary balancing, we will obtain negative values only.

The Balancing Procedure
An important step in applying the functional expansion method is related, as in almost all expansion methods, to the balancing procedures, which allow the limit of the expansions to be found. As we already mentioned, the functional expansion method supposes two different expansions: one in terms of various powers of the derivatives G (ξ), and one in the chosen form of the functionals P i (G). This fact automatically leads to two balancing procedures: the first one allows the maximal value of m in (12) to be found, while the second one leads to the possible forms of the functionals P i (G) and allows their degrees to be determined, n(P i ), as defined by (16). When n(P i ) are known, n(N i ) and n(D i ) can be fixed, when the representation in (15) is considered.
Speaking about balancing, an important issue to also be considered is the mathematical form of the equation to be solved. To illustrate how effectively all these balancing procedures are working, let us consider that the ODE (3) to be solved has the following general form: This class of NODEs leads to many important nonlinear 2D equations. One of the examples, which was intensively tackled in this work, is the Dodd-Bullogh-Mikhailov equation. In this case, the attached NODE has the following form: Other important equations lead to NODEs corresponding to B(u) = 0. Such equations are as follows: • the KdV equation • the cubic nonlinear Schrödinger equation • the nonlinear Klein-Gordon equation • the ZKBBM equation The previously mentioned equations could have more general nonlinear terms. For example, in the Klein-Gordon equation considered in quantum field theory, other nonlinearities can appear [37][38][39][40].
To keep the discussions as general as possible, we consider that in (17) A(u), B(u) and E(u) are polynomials of the following form: We start by discussing the last aspect mentioned: a balancing procedure imposed by the form of the equation to be solved. We suppose the functions in (23) are known, that is the limit of each expansion, n A , n B and n E , are known. By introducing (12) in (17) a system of nonlinear ordinary differential equations for the functionals {P i (G), i = 0, 1, · · · , m} is generated by equating the coefficients of the terms with the same powers in G to zero. This is called the determining system and is used for finding the degrees attached to the functionals, following (16).
In the first step, we see what values could have the summation limit m appearing in (12), as function of the parameters n A , n B and n E . They can be obtained through a first balancing procedure among the term with the highest derivative and the terms with the biggest order of nonlinearity. The maximal order of derivation in G is generated by the first term from (17). It has the following form: Depending on the degrees of the polynomials B(u) and E(u), the terms that generate the highest nonlinearity are as follows: If, for example, n E > n B , the first balancing has to be made between (24) and (26). It allows us to find the summation limit m in (12), leading to m(n A + 1) + 2 = mn E : Imposing m ∈ N, we conclude that m can take two integer values only: either m = 2 (then n E − n A = 2) or m = 1 (for n E − n A = 3). The case n E − n A = 1 asks for special consideration, which is not the object of analysis in the current paper. Such a situation appears in the Chafee-Infante model: Let us go further by trying to get the conditions the functionals {P i , i ≤ m} have to satisfy, for an already fixed m. For this purpose, the equation in P m , the functional of the maximum degree with m given by (27), contains only two terms and it has the following form:P We introduced the notation α = e n E a n A , the ratio of the coefficients appearing in front of the maximal order terms in the expansions from (23). This is a constant with a known value for a given model. Let us take P m (G) as a rational polynomial expression, having the form (15). The compatibility requirement for (28) allows, by applying a second balancing procedure, the degree n(P m ), attached to P m by (16), to be effectively determined. More precisely, it leads to a relation among the already known m and the degree n(P m ): Considering the case described by (27) with m = 2, which will be tackled below, as m is positive, we deduce that n(P m ) has to be negative.
Until now, m and the functional P m have been determined. The other functionals {P i , i < m} appearing in the solution (12) can be determined considering the other equations of the determining system that were generated when (12) was introduced in (17). Similar reasonings as before lead to negative values for all the functionals' degrees: This will be proven in the next section, using one of the equations identified as belonging to (17) as an example. (15) being quite similar to that used in [27], the relation (29) leads our approach to n(N ) <n(D), while [27] shows the opposite: the numerators have higher degrees than the denominators.

Remark 2.
It is important to note that (29) and (30) fix only the difference between the two degrees n(N i ) and n(D i ). It is clear that there are many solutions that can be considered and that there are many possible choices of the type in (15) for the same functional P i (G). For m = 1, (30) imposes n(P i ) = −1, which can be achieved considering, for example, n(N i ) =1 and n(D i ) =2, but also n(N i ) =2 and n(D i ) =3.

The Example of the Dodd-Boullogh-Mikhailov Equation
To see explicitly how the previous assertions functioned we considered a specific model of 2D NPDEs, leading, when the wave variable (2) was introduced, to a particular case of (23), namely to (18); this is the Dodd-Bullough-Mikhailov (DBM) equation, with the following form: With the change in variable u(x, t) = e w(x,t) the previous equation takes the following form: This is an important equation, with many applications in hydrodynamics and quantum field theory. Various types of periodic, hyperbolic or rational solutions, of traveling wave or of soliton types, were pointed out, using methods such as the tanh method [41], the exp-function method [42] or the G /G method [43,44]. Here the equation is investigated using the functional expansion method and, as we will see, this approach allows the recovery of all the mentioned solutions, and, moreover, enables new solutions, larger that the G /G solutions, for example, to be found.
As we already mentioned, the first step of the functional expansion method consists of the reduction of (32) to an ODE using the wave variable (2). This reduction leads to (18), that can be re-written as follows:

The Determining System for the Functionals
Let us look for solutions of the DBM equation of the form (12). We use the expression in (33) and we apply the first balancing procedure for determining the number of terms to be considered in the expansion. Taking into account the term with a second order derivative and third order nonlinearity, the balancing leads to m = 2, so the sought DBM solution has the following form: Finding the DBM solutions allows P 0 , P 1 , P 2 to be found as functionals depending on the solutions G(ξ) of an auxiliary equation. We choose here, as an auxiliary equation, a general second order differential equation: From (33)- (35), equating the coefficients of the different powers of G to zero, we obtain by hand, but also using Wolfram Mathematica, a determining system of seven ordinary differential equations for P 0 , P 1 , P 2 : −V ..
−µV (42) from the previous generating system can be rewritten as follows:

Remark 3. The last Equation
As shown below, this equation leads, in all the cases, to a constraint showing that the parameters λ and µ from the auxiliary equation cannot take any values; they are related to each other and to the wave velocity V. The same constraint also appears when G /G solutions are considered.

Remark 4.
A first attempt at directly solving the previous obtained system would probably lead to the most general solution accepted by the DBM model. It is quite easy to verify, for example, that the first Equation (36) accepts the following as a solution: For C 1 = V and C 2 = 0 the equation becomes Using (45), (37) can be solved, obtaining the solution for P 1 : Unfortunately, this approach of finding DBM solutions by directly solving the determining system (36)-(42) fails when trying to find P 0 from the remaining Equations (38)- (42). It seems that it is not possible, at least for the DBM model, to obtain general P 0 , P 1 and P 2 that are compatible with the whole system. This is why another approach is needed to find solutions, with the functionals chosen as in (15).
An important step is determining the limits n(N i ) and n(D i ) in the expansions of the numerator and the denominator of each {P i , i = 0, 1, 2}. For this purpose the second balancing procedure is used, applied this time to the determining Equations (36)- (42). Taking into account that for DBM we obtained m = 2, from (29) we obtain n(P 2 ) = −m = −2. Similarly, we conclude that the degrees attached to the functionals {P i , i = 0, 1, 2} by (16) have to be as follows: n(P 2 ) ≡ n 2 ≡ n(N 2 ) − n(D 2 ) = −2.
The limited constraints (47) offer a large freedom in the choice of the mathematical form of the functionals P 0 , P 1 and P 2 . Correspondingly, larger classes of DBM solutions may be generated through the functional expansion method. As a proof, in the next subsection we analyze three choices that are more general than those considered in the G /G method.

Remark 5.
Let us mention again that, in all the three cases, the simplest choice corresponds to the following: The DBM solution (34) has exactly the form given by the G /G-approaches: It is obvious that the choices (48)-(50) are more general.

Examples of DBM Solutions Generated through the Functional Expansion
We will now show how our proposed method effectively functioned, to see if more general solutions, such as those arising in the G /G method, could be generated. The procedure was quite simple and obvious: we introduced the chosen forms for the functionals P i in the determining system (36)- (42), taking into consideration the explicit form of the auxiliary Equation (35). A set of algebraic equations arose, relating the parameters { ω ij , π ij } from P i with the wave velocity V from the main Equation (33) and with the parameters {λ, µ, ρ} from the auxiliary equation. All the compatible solutions of this algebraic system led to solutions for the functionals {P i , i = 0, 1, 2}, and, implicitly, for the DBM Equation.
We have to keep in mind that the solution of (35) could be written as in [45]: Depending on the relation between λ and µ, we have three different situations: (i) If λ 2 − 4µ < 0 we have the following: Here, as well as in the forthcoming expression, we used the notations C 1 = Ae −iϕ and C 2 = Ae iϕ , respectively. (ii) If λ 2 − 4µ > 0, the solution (53) could be written as (iii) If λ 2 − 4µ = 0 the solution (53) is:

Examples of Solutions in Case I
The functionals {P i , i = 0, 1, 2} have numerators of degree zero: n(N 0 ) = n(N 1 ) = n(N 2 ) = 0. To observe the constraints (47), the denominators of the functionals have to be n(D 0 ) = 0; n(D 1 ) = 1; n(D 2 ) = 2. Choosing simplified notations for the coefficients appearing in (48)-(50), we consider the following: We note that with the choice in (57), Equation (42) in the determining system leads to the following constraint: Many DBM solutions can be generated with these choices. Some of them correspond to the already-reported solutions, obtained through the G /G-method. For example, one of the solutions accepted by the determining system (36)-(42) is of the following form: This corresponds to the case ρ = 0, λ = 0, and it leads to the solution of (33) of the following form: On the other hand, even observing the constraint in (58), non-standard solutions of the determining system appear as follows, for example: The corresponding solution of (33) becomes the following: Other solutions, apparently more complex than (62), are as follows: Although, when (58) is imposed, these solutions take the form of (61).

Examples of Solutions in the Case II
When the functionals {P i , i = 0, 1, 2} have numerators of degree one, n(N 0 ) = n(N 1 ) = n(N 2 ) = 1, the constraint in (47) asks for n(D 0 ) = 1; n(D 1 ) = 2; n(D 3 ) = 3. Again, considering simplified notation for the coefficients in (48)-(50), we choose the following: We also consider here that ρ = 0. The algebraic equations generated by the determining system (36)- (42) lead to the relations among the parameters {a, b, c, d, e, f , h, j, n, p, r, V, λ, µ} and allow the functionals P 0 , P 1 and P 2 to be in a simpler form: Compatibility conditions impose, in this case too, the following supplementary constraint: Similarly to (58), it restricts the possible values of the parameters in the auxiliary equation. As the velocity is a real quantity, and, as we have seen, the solutions of the auxiliary equation ask λ 2 − 4µ to also be real (positive, negative or zero), we retain from (68) the following situation: We note that only two situations, (54) and (55), can be fulfilled, considering adequate values, positive or negative, for the velocity V. Correspondingly, we have to consider only these two types of solutions, harmonic and hyperbolic, for the auxiliary equation. There are no realistic velocities leading to λ 2 − 4µ = 0.
Let us also note that for ρ = 0, the functionals (65)-(67) become the following: The final DBM solution of Equation (33) becomes, in this case, the already-known expressions provided by the G /G approach: Another remark is that, when (68) is observed, the expressions (66) and (67) take the simple forms from (61) mentioned in Case I.

Examples of Solutions in Case III
Consider now the case when all the functionals P i have identical quadratic denominators: n(D 0 ) = n(D 1 ) = n(D 2 ) = 2. The relation in (47) imposes that n(N 0 ) = 2; n(N 1 ) = 1; n(N 2 ) = 0. We may choose P 0 , P 1 , P 2 as having the following forms: The procedure mentioned before leads, in this case too, to non-standard G /G-solutions. It is interesting that, again, Equation (42) is fulfilled if and only if the wave velocity is related to λ and µ from (35) by a relation of the same form as in the two previous cases: In fact (69) and Equation (72) now also take the following simplified forms: Again, because of (69), we should consider only the harmonic and the hyperbolic solutions of auxiliary Equation (35), that is (54), respectively (55). For example, for negative velocities, we have λ 2 − 4µ > 0, and for positive velocities, λ 2 − 4µ < 0.
Below some comments on DBM solutions are giventhat were obtained through the functional expansion method, in the three examples that were considered before; there were some similarities that did not depend on the chosen form of the functionals P 0 , P 1 , P 2 . Equation (42) generates, in all the cases, the following constraint: In all the cases the general solutions obtained could be, at the end, reduced to the same expressions: Using the expressions in (34), they led to DBM solutions that were different from those that the G /G approach generated.
Apparently, the non-standard G /G-solutions (75) are related to the factor ρ = 0 in (35). It is important to note that such non-standard solutions appear even if ρ = 0. It is simple to check, for example, the following two solutions for the determining system (36)- (42) corresponding to ρ = 0: (76) These are also different from what was obtained by the G /G approach.

Recovering the Main Types of DBM Solutions
The DBM solutions obtained using functional expansion had the general form (34). We effectively wrote down a few of these solutions, using the expressions (75)-(77) for the functionals {P i , i = 0, 1, 2}. It was interesting to note that, whatever expression was used, we obtained quite similar DBM solutions. They depended on the wave velocities and all the parameters appearing in (75)-(77) were captured in two other parameters that are denoted below by C 1 and C 2 , respectively. As already mentioned, only two cases arose and they corresponded to the solutions of (55) and (54), respectively, of auxiliary Equation (35). We prove here that these two cases, corresponding to negative and, respectively, positive wave velocities, practically allowed the recovery of all the important types of DBM solutions.
For negative velocities, V = − 3 λ 2 −4µ < 0, the auxiliary equation admits the solution of (55), and it leads to the DBM solution (C 1 and C 2 are integration constants): In Figure 1a this solution is represented for V = −1 and for any C 1 = C 2 . It has the form of a bright soliton, which is the type of solution already reported in the literature for the DBM model. Its specific mathematical form is as follows: For V = −1 and any C 1 = −C 2 , the solution is plotted in Figure 1b. It is the typical dark soliton accepted by DBM and can be re-written as follows: If we considered bigger values for the two constants, C 1 and C 2 , the solution profile (78) changed. It took the form of periodic peaks propagating in time and along the x-axis. The peaks could have unbounded amplitudes and a periodicity depending on the effective values of C 1 and C 2 . This behaviour is illustrated by the two specific solutions plotted in Figure 2. In principle, bigger values of the constants C 1 , C 2 led to decreased wave periods and amplitudes. We noticed that the amplitudes changed from ∼10 30 in Figure 1b to ∼10 4 in Figure 2. For positive velocities, V = − 3 λ 2 −4µ > 0, the auxiliary equation admitted the solution in (54). In this case the DBM solution took the following form: Considering V = 1, the solution (81) for C 1 = 0, C 2 = 1 is plotted in Figure 3a. It had the form of many propagating periodic waves. For bigger values of the two constants, C 1 = 5, C 2 = 7, the wave amplitudes and frequencies decreased, as can be seen in Figure 3. For C 1 and C 2 with opposite signs, the solution had a similar shape: many waves propagating along the axis. From Figure 4a,b, made for different values of C 1 , C 2 , we note that there was not a high dependency on the values of these constants. This was quite normal considering the mathematical form of the solution.

Conclusions
This paper presented in detail how the functional expansion method proposed in [32] functioned for a large and important class of equations that can be expressed as in (17). Such nonlinear equations have important applications in various fields, such as optics and plasma physics [46][47][48], for example. Our claim was that this approach for solving NPDEs was more general than almost all the others based on the use of an auxiliary equation. Two such approaches were specially considered: the Kudryashov method, which is suitable when first order auxiliary equations are considered and the G /G method, which is the traditional approach in the case when the focus is on second order auxiliary equations. Compared with previously published papers, the novelty this paper brings was related to the explicit presentation of the balancing procedure that, in the case of the functional expansion method, required a double balance: the first one gave the maximal term in the expansion of (12), and the second one was used for determining the functionals P 0 , P 1 , P 2 , as is explained in Section 3.1. The application of the functional expansion to the DBM model represented another novelty of this paper.
The method was based on expansions of the type in (10), or, more exactly, of the type in (12). These were in fact the most general possible forms of solutions and they included almost all the choices used in various approaches to the direct finding of exact solutions for nonlinear differential equations. The method presented many advantages, one of them being that it generalized other approaches to the direct solving of NPDEs. Here we considered the Kudryashov and (G /G) methods [28]. The choice of (15) is similar to what Kudryashov method considers. Practically, the functional expansion approach extended the Kudryashov approach to second order auxiliary equations, and it allowed more general solutions, as in the (G /G) approach, to be obtained. This was another important merit of out method and it was illustrated for the DBM equation using the non-standard solutions of type (75)-(77). Expressions containing the G /G ratio appear now in the most natural way, as particular sub-cases of more general solutions. It was true that the non-standard form of solutions were limited to first order denominators; this introduced a limitation to our method, at least for the DBM model.
Another important issue approached in the paper was related to the balancing procedures that were traditionally applied to limit the number of terms considered in the expansions. It was pointed out that the functional expansion asks for two different balancing procedures: one following the powers of G and a second one following the powers of G. The connection between the two, as well as the relation with the form of the equations, were investigated for equations belonging to the class in (17). The outcome expressed through (27) was quite important for investigating equations such as KdV, nonlinear Schrödinger, Klein-Gordon, KMN [33] or Benjamin-Bona-Mahony. Limits of the method in investigating special types of equations, such as Chaffe-Infante or Fisher, for example, were also mentioned. How these limitations could be overcome using alternative approaches will be tackled in future works. Acknowledgments: Three of the authors (C.N.B., R.C. and R.E.) acknowledge the mobility grants received from the H2020 Project "Dynamics", 2017-RISE-777911, as well as the support they obtained in the frame of the NT-03 Agreement between SEENET-MTP and ICTP.

Conflicts of Interest:
The authors declare no conflict of interests.

Abbreviations
The following abbreviations are used in this manuscript: NPDEs Nonlinear partial differential equations NODE Nonlinear ordinary differential equations ODE Ordinary differential equations KdV Korteweg de Vries KMN Kundu-Mukherjee-Naskar DBM Dodd-Bullough-Mikhailov