On Solving the Nonlinear Falkner–Skan Boundary-Value Problem: A Review

: This article is a review of ongoing research on analytical, numerical, and mixed methods for the solution of the third-order nonlinear Falkner–Skan boundary-value problem, which models the non-dimensional velocity distribution in the laminar boundary layer.


Introduction
The Falkner-Skan equation, a third-order ordinary differential equation over a semiinfinite domain, was originally developed in 1931 by Falkner and Skan [1] from a similarity analysis on the steady, two-dimensional, boundary-layer flow of an incompressible, viscous fluid over a wedge [2]. As it turns out, the solution of the Falkner-Skan equation represents the laminar boundary layer over an infinite wedge whose vertex angle is βπ for 0 ≤ β ≤ 2, as illustrated in Figure 1. The case β = 0 corresponds to the flow past a flat plate and leads to the well-known Blasius equation [3], and the case β = 1 corresponds to what is known as the stagnation-point flow or the Hiemenz flow [4]. The similarity analysis begins with the continuity and the momentum equations that are described by ∂u ∂x subject to the boundary conditions u(x, 0) = v(x, 0) = 0, and u(x, y) → U(x) as y → ∞.
In (1)- (3), u and v are the velocity components in the x and y directions of the fluid flow; µ is the viscosity; and U(x) is the velocity at the end of the boundary layer. A stream function Ψ(x, y) is defined so that the continuity Equation (1) is automatically satisfied. Accordingly, it must be true that u = ∂Ψ ∂y , and v = − ∂Ψ ∂x .
derived from the corresponding transformation of (3). In (5) and (6), the primes denote differentiation with respect to η. The nonlinear, third-order, ordinary differential Equation (5) is known as the Falkner-Skan equation. Together, (5) and (6) will be referred to as the Falkner-Skan boundary-value problem. Solutions of (5) and (6) are characterized by the numerical value α = f (0). Many industrial processes such as the drawing of plastic films, metal spinning, cooling of a metallic plate, and aerodynamic extrusion of plastic sheets require the solution of the Falkner-Skan equation and several of its variants applied to different geometries and corresponding boundary conditions [5]. Thus, there has been a considerable research interest in solving the Falkner-Skan equation during the past few decades. Analytical investigations of this problem have been aimed at providing existence and uniqueness results, finding exact, nearly-exact, or approximate analytical solutions. Computational approaches include a spectrum of methods ranging from the traditional finite-difference and finite-element methods to applying neural networks. The purpose of this review article is to provide a succinct summary of these various approaches that have been used for solving this important nonlinear boundary-value problem.

Materials and Methods
For the purpose of this review, we will be specifically investigating the problem as posed by (5) and (6). As this is a review article focused on the existing literature, this section will be devoted to describing the methods and/or approaches used by previous researchers. Accordingly, this section will include two main subsections, the first one reviewing analytical treatments, and the second exploring the computational approaches documented thus far in the literature. The discussion of analytical treatments will include existence and uniqueness results, exact or nearly-exact analytical solutions, and approximate analytical solutions. The discussion of computational approaches will include the description of the different classes of methods currently available.

Analytical Treatments
The first result on the existence and uniqueness of a solution to (5) and (6), with the additional assumption of f > 0, is due to Weyl [6] and this result was proved using a fixed-point theorem. Later, Coppel [7] assumed 0 < f < 1, considered three-dimensional trajectories in the phase-space, and established the existence and uniqueness result of Weyl [6], along with many more details on the asymptotic behavior of the solutions with β ≥ 0 for various scenarios related to β in (5). In addition, Coppel [7] also showed that α = f (0) is an increasing function of β. Tam [8] removed Coppel's assumption of 0 < f < 1 and established the existence and uniqueness result using a shooting argument. Craven et al. [9] proved the existence and uniqueness of solutions for 0 ≤ β < 1 by arguing that in this range, the condition 0 < f < 1 is satisfied. Solutions of (5) and (6) which have f (η) < 0 for η > 0 represent what are known as reverse flows. Stewartson [10] was the first to obtain reverse flows for β ∈ (−0.1988, 0). The existence and uniqueness of the reverse flows obtained by Stewartson was proved by Hastings [11]. In related work, Craven et al. [12] carried out numerical experiments that suggested the existence of reverse flows for β > 1. Other works on existence and uniqueness of solutions to (5) and (6) may also be found in Rosenhead et al. [13] and Hartman [14]. Veldman and Van der Vooren [15] have extended Coppel's work for β < 0. The work of Oskam et al. [16] is also focused on the study of (5) for β < 0. In place of the boundary conditions in (6), Oskam et al. [16] consider two scenarios, one corresponding to f (∞) = 1, and one corresponding to f (∞) = −1, which are closely related as β → 0 − . Moreover, for β < −1, these two scenarios allow multiple solutions which are characterized and distinguished using the number of zeros of f − 1. Additionally, for β < −1, this work provides a periodic solution as well, and a more general treatment of periodic solutions is presented by Botta et al. [17]. Rajagopal et al. [18] have established non-similar solutions for the boundary layer flow of a homogeneous incompressible fluid of second grade past a wedge placed symmetrically with respect to the flow direction and discussed the variation of the skin-friction with respect to the non-Newtonian parameters. In particular, by expressing the solution as a power series in terms of a parameter ε that scales velocities with the square root of the Reynolds number, they obtain the zeroth-order and first-order solutions. Their zeroth-order solution coincides with the classical Falkner-Skan flow which we have explored in this review. They have numerically computed the second-order solution and analyzed it more in-depth as it relates to variations in the parameter ε. Massoudi et al. [19] have investigated the effects of the non-Newtonian nature of a fluid on the shear stress at the wall as well as the velocity profiles for different injection and suction rates and different values of the wedge angle β. Additional existence and uniqueness results for the solution of the Falkner-Skan equations are also given by Padé [20].
More generally, with (5) written in the form subject to the boundary conditions (6), several instances corresponding to specific choices of values for β 0 and β have been explored analytically.

The Pohlhausen Flow
For (β 0 , β) = (0, 1), (7) leads to what is known as the Pohlhausen flow [21]. This instance, described by has the exact solution given by It is sometimes used to validate numerical methods [22].

The Blasius Flow
With (β 0 , β) = ( 1 2 , 0) (7) takes the form and leads to the well-known Blasius probem [3]. While a closed-form solution is not available for (10), it has been subject to extensive analytical treatment as well as numerical investigation. We briefly describe some selected analytical treatments here and postpone the discussion of the numerical investigations to a later section. For the case of (10), Blasius [3] proposed the power series solution with and α = f (0). This solution requires knowledge of α = f (0) and, several works [23,24] have pointed out that the above series is valid only for small values of η. As the search for an approximate analytical solution of (10) continued, Ahmad and Al-Barakati [25] modified the [4/3] Padé approximant for the derivative to satisfy its asymptotic behavior and obtained the solution Beginning with an approximate solution of the form f (η) = 1/(Aη + B) that contains the undetermined coefficients A and B, Bougoffa and Wazwaz [26] use the differential Equation (10) and the boundary conditions (6) to obtain the approximate analytical solution Other works that provide asymptotic approximants to the solution of (5) and (6) include Belden et al. [27].

The Crocco-Wang Transformation
Another interesting line of investigation of the Falkner-Skan and Blasius flows has evolved due to the suggestion of the transformation that originally appeared in Crocco [28], and later independently used by Wang [29]. Under the transformation (15), the third-order boundary-value problem (10) along with the boundary conditions (6) on the semi-infinite domain [0, ∞) takes the form which is a second-order boundary-value problem on a finite interval. Note that in (16) the primes denote differentiation with respect to the variable x. In this form, the quantity y(x = 0) corresponds to α = f (η = 0) that characterizes the solution of the Blasius problem. Brighi et al. [30] provides a historical review of the research on the Blasius problem. Several analytical approaches to analyzing and solving (16) are currently available.
The work of Varin [31], where y(x = 0) is referred to as the Blasius constant, examines the Crocco Equation (16) and provides the following results: (1) The non-existence of singularities other than the x = 1 in (16) (which represents η → ∞ in (6)); (2) an explicit expression for the Blasius constant as a sum of a series of rational numbers; (3) a recursive formula for the coefficients of the series; and (4) a numerical procedure for computing the Blasius constant to high accuracy.
Under the Crocco-Wang transformation (15), the Falkner-Skan problem (5) and (6) becomes , and y(1) = 0, (17) in which the primes denote differentiation with respect to the variable x. Note that y(x = 0) corresponds to α = f (η = 0) also for the transformed Falkner-Skan equation as it did for the Blasius equation. The existence of a unique positive solution to a variant of the instance of (17) corresponding to the Homann [32] flow, or (β 0 , β) = (1, 1 2 ), has been demonstrated in Shin [33].

Other More General Analytical Approaches
Adomian Decomposition Method (ADM). Rewriting (5) . First, f 0 (η) in the infinite series expansion is obtained as the solution of the homogeneous equation L f 0 = 0, which yields In order to obtain f n (η), the ADM defines the Adomian polynomials as and sets According to Alizadeh et al. [34], the ADM is more realistic compared to other methods that may use linearization in some form or another, or assume some weak nonlinearity, or simplify the physical problem. On the other hand, the ADM involves a considerable amount of symbolic manipulation in order to obtain more terms in the infinite series expansion. For the Falkner-Skan problem, the results reported in [34] indicate that ADM required from 2 to 30 terms in the infinite series expansion as the wedge angle was varied from 0 • (flat plate) to 180 • (stagnation flow).
Homotopy Methods. A variety of methods, known by the names of Homotopy Asymptotic Method (HAM) [35], Optimal Homotopy Asymptotic Method (OAHM) [36][37][38], and Homotopy Perturbation Method (HPM) [39,40] have also been used to obtain approximate analytical solutions to the Falkner-Skan problem (5) and (6). These methods express the problem in the form where L is an easily invertible linear part of the operator representing the Falkner-Skan problem, N the entire nonlinear operator that represents the Falkner-Skan equation, and p a parameter that varies from 0 to 1. The form (20) is set up such that the solution corresponding to p = 0 can be easily obtained by inverting the linear operator, and the solution corresponding to p = 1 is the desired solution of the Falkner-Skan problem, both satisfying the boundary conditions (6). In order to advance the solution from p = 0 to p = 1, the desired solution is assumed to be of the form Beginning with f 0 (η) obtained as the solution of Lφ(η) = 0 subject to the boundary conditions (6), plugging the form (21) into (20), differentiating with respect to p as many times as desired and setting p = 0 yields a succession of equations to solve for f k (η).
An efficient implementation of HAM by the name Spectral HAM (SHAM) [41] has also been proposed, and this method obtains the initial solution using collocation at Chebyshev points after suitably transforming the semi-infinite domain into the interval [−1, 1]. It is customary to call p the embedding parameter. In addition to p, the HAM approach allows for the use of another parameter h known as the convergence control parameter.
As evident from the above discussion, most studies focused on the analytical treatment of (5) and (6) or its variants have demonstrated a solution's existence and uniqueness, developed semi-analytic, asymptotic solutions for special cases, or discovered ranges of validity for the boundary-layer parameters and similarity variables [5,6,8,13,14,20,42]. However, a general closed-form solution is still unavailable for this two-point boundaryvalue problem in spite of the efforts dedicated to its investigation. Therefore, approaching its solution numerically has been the most common and valuable route.

Numerical Treatments
Several different numerical approaches have been used to investigating the Falkner -Skan problem (5) and (6), and they are in the categories of shooting, finite-differences, collocation, and other methods. The first computational treatment reported in the literature is due to Hartree [43]. While the details of the numerical approach are not provided in this work, Hartree [43] showed that physically relevant solutions of (5) and (6) exist only in the range −0.198838 ≤ β ≤ 4/3. Smith [44], Cebeci and Keller [45], and Na [46] have used shooting and invariant imbedding to solve the Falkner-Skan problem directly on the semi-infinite domain. The methods that came about after these early treatments mainly focused on making the computational methods more efficient and/or skillfully dealing with the boundary condition at infinity. Various special cases of (5) and (6) have been used to validate the effectiveness of the methods developed.
Coordinate transformations and appropriate change of variables have been used for the purpose of mapping the physical domain [0, ∞) to a computational domain like [0, 1] or [−1, 1]. It is also common practice to truncate the physical domain where η → ∞ by a finite value η = η ∞ which is determined as part of the solution, and impose the boundary condition as η → ∞ at the truncated boundary η = η ∞ . Once transformed to a suitable computational form, methods vary depending on the specific way in which numerical computations are carried out. Among the different techniques used for this purpose are, shooting; finite-differences; and collocation methods. Along with these fundamental differences in the approach, the methods also vary based on the specific techniques used for the computational problems that result. In this section, we will review these approaches and the associated computational techniques.

Shooting
The term shooting in numerical analysis refers to the idea of solving a two-point boundary-value problem iteratively by repeatedly solving an initial-value problem aimed at obtaining the value such that the boundary condition as η → ∞ is satisfied. With an initial guess for α, the shooting process involves marching the solution from η = 0 to η = ∞, and adjusting or improving the guess for the value of α appropriately so that the boundary condition at η = ∞ is satisfied. There are several ways in which the methods that employ the shooting approach may differ. These include the initial guess for α, the choice of the method that advances the solution from η = 0 to η = ∞, and the manner in which the guess for α is adjusted. As mentioned earlier, the semi-infinite physical domain [0, ∞) is first mapped to a finite computational domain like [0, 1] using an appropriate coordinate transformation. For this purpose, a finite boundary η = η ∞ is first introduced and the value of η ∞ will be determined as part of the solution. It is customary to introduce the "asymptotic condition" d 2 f dη 2 = 0 as η → ∞ to help determine the additional unknown η ∞ . Thus, with the truncated domain [0, η ∞ ], the conditions will be imposed. Next, we observe that under the coordinate transformation ξ = η/η ∞ [47], the problem (5) and (6) becomes subject to the boundary conditions At this point, applying the the change of variables transforms the single third-order differential equation to (24) to the system of three firstorder differential equations while the boundary conditions (25), (26) and (22) become Additionally, the conditions (23) now take the form u = 1 and v = 0 at ξ = 1.
Note that (29) and (30) together represent an initial value problem. With the correct value of α, the initial value problem (29) and (30) should produce u and v that satisfy the conditions (31). Conversely, in order to obtain the correct value for α, the initial value problem (29) and (30) is first solved with an initial guess for α in (29), the values of the u and v at ξ = 1 based on this solution are used to improve upon the guess for α. This process is repeated until the conditions (31) are satisfied. The following is a listing of the works that have used the shooting approach along with a description of how they differed from the others using the same approach: 1.
The work of Cebeci and Keller [45] is the first to employ the shooting approach to solve (5) and (6). They work directly with the physical domain η ∈ [0, ∞), but march their solution from η = 0 to η = η ∞ , where they assume η ∞ is "sufficiently large." They use the Newton's method to adjust the initial condition.

2.
Cebeci and Keller [45] also employ the idea of parallel shooting to obtain the solution. This approach involves dividing the domain into several subintervals and applying shooting in each subinterval. The motivation for this approach comes from the desire to avoid any instabilities that may be caused while marching the solution using an initial-value method. The authors report successfully obtaining three significant digits in two iterations when shooting was carried out over three subintervals. 3.
Asaithambi [47] uses the classical, fourth-order Runge-Kutta method to march the solution from ξ = 0 to ξ = 1 and the secant method to iteratively obtain improved values for α. Two initial approximations for α were needed. 4.
Salama [48] uses a single-step method of order 5 to march the solution from ξ = 0 to ξ = 1, coupled with a third-order rootfinding method called the Halley's method [49] to compute α and the secant method to compute η ∞ . It is also important to note that the solutions of (5) and (6) are obtained for β values as large as 40.

5.
In another work, Asaithambi [22] uses recursive evaluation of Taylor coefficients to march the solution from ξ = 0 to ξ = 1 and the secant method to iteratively obtain improved values for α. Two initial approximations for α were needed, and higher orders of accuracy, as high as 7, were achieved with ease. Moreover, this work explores (5) and (6) for β values as large as 40 as well and compares the results obtained previously by Salama [48]. 6.
Zhang and Chen [50] used the Runge-Kutta (4,5) formula to march the solution from ξ = 0 to ξ = 1, and the Newton's method to improve the guess for α, and hence required only one initial guess to get started. 7.
Sher and Yakhot [51] shoot from the end corresponding to η → ∞ instead of from η = 0. They choose an arbitrary value for f (η) as η → ∞, a value close to 1 for f (η) as η → ∞ and estimate f (η) as η → ∞ using a simple analysis of the asymptotic behavior of the solution.

Finite Differences
The finite-difference approach is often used as a more efficient alternative to the shooting approach, as the use of finite differences does not warrant the repeated solution of initial-value problems as does the shooting approach. The following is a listing of the works that have used finite differences, along with a description of how they differed from the others using the same approach:

1.
Asaithambi [52] truncates the infinite domain at η = η ∞ , uses the transformation ξ = η/η ∞ as done previously, and sets u = 1 η ∞ d f dξ , to transform the problem (5) and (6) on the physical domain η ∈ [0, ∞) to the system on the computational domain ξ ∈ [0, 1], subject to the conditions In order to discretize (32)- (35), the interval ξ ∈ [0, 1] is subdivided into N subintervals [ξ j , ξ j+1 ] for j = 0, 1, · · · , N − 1 at ξ j = jh with h = 1/N. Letting f j and u j denote the approximations to f (ξ j ) and u(ξ j ) respectively, Equation (32) is discretized as and Equation (33) is discretized as The boundary conditions (34) translate to and the conditions (35) translate to Newton's method is used for the iterative solution of the nonlinear system (37), Equation (36) is used to update the values of f j , and (39) is used to update η ∞ . The solution of the nonlinear system by Newton's method involves the repeated solution of a linear system with a tridiagonal coefficient matrix. In related work, Asaithambi [53] uses the transformation ξ = e −η to map the physical domain [0, ∞) to the computational domain [0, 1] and discretizes the transformed equations and boundary conditions using second-order finite differences. This approach does not truncate the physical domain but still requires the iterative solution of a nonlinear system that involves the solution of a tridiagonal linear system during each iteration.

2.
Salama and Mansour [54] develop an unconditionally stable, fourth order finite difference method that involves only four grid points and apply it to the nonlinear Falkner-Skan problem (5) and (6). Their development uses the method of undetermined coefficients involving the unknown function values at points η j , η j−1 , η j−2 , η j+1 to derive the discretization for the third derivative by setting the local truncation errors of orders 0, 1, 2, and 3 to zero. The details of the derivation are too technical to be included here and the reader is referred to the original publication [54]. In related work, Salama and Mansour [55] develop a finite difference method of order 6 for solving steady and unsteady two-dimensional laminar boundary-layer equations. This method is also unconditionally stable. The authors illustrate the effectiveness of this method by solving the unsteady separated stagnation point flow, the Falkner-Skan equation and Blasius equation.

3.
Duque-Daza et al. [56] express (24) as and use N − 1 subintervals with h = 1/(N − 1) and ξ j = (j − 1)h for j = 1, 2, · · · , N − 1. With , and f (ξ i ) respectively, the authors use the fourth-order difference formulas for the first and second derivatives of f and the third-order formula for the third derivative. Plugging (41)-(43) into (40) yields a nonlinear system of N − 5 algebraic equations in the N − 1 unknowns f 1 , f 2 , · · · , f N−1 . The boundary conditions are discretized using Next, plugging f = η ∞ and f = 0 into (40) for ξ = 1 yields which is discretized as Then the Equation (40) itself is evaluated at ξ = 0 to obtain f 1 + βη 3 ∞ = 0, which is discretized as Finally, the condition f N = 0 is used to solve for η ∞ . For this purpose, the secondorder derivative f N is discretized as and the secant method is used to obtain updated values of η ∞ from one iteration to the next. In the same article [56], the authors also report the use of a higher order compact finite difference method to solve (5) and (6).

Collocation Methods
The collocation approach for solving a boundary-value problem consists of selecting a solution from a finite dimensional space of candidate solutions by requiring that such a solution satisfy the governing differential equation(s) at a select set of points in the domain called collocation points. More specifically, the desired solution from the finite dimensional space is expressed as a linear combination of a set of basis functions with undetermined coefficients. Requiring this solution to satisfy the differential equation and the boundary conditions together reduces the boundary-value problem to the problem of solving a system of algebraic equations in the undetermined coefficients. While most collocation methods for solving the Falkner-Skan problem (5) and (6) attempt to obtain polynomial forms for the solution, they differ based on the choice of basis functions and the techniques they employ in order to translate the given boundary-value problem to the system of algebraic equations. We review a few such methods here.

1.
Parand et al. [57] consider (5) subject to the boundary conditions in which f (0) = γ is different from the corresponding condition in (6). The authors use rational Legendre functions as their basis functions, and the correspondingly transformed Gauss-Legendre points as the collocation points. It is important to emphasize here that their approach does not require the problem to be transformed to a suitable form on a finite domain. The rational Legendre functions {R n (η)} are defined by in which {P n (y)} are the well-known Legendre polynomials. While P n (y) form a family of polynomials orthogonal on the interval [−1, 1] with respect to the weight function σ(y) ≡ 1, R n (η) form a family of rational functions that are orthogonal on [0, ∞) with respect to the weight function Moreover, in (52) and (53), L is a constant known as the map parameter. While setting L = 1 is most common, guidelines for the choice of value for L are given for rational Chebyshev functions in [58] that apply also to the Legendre family. The rational Legendre functions satisfy the relation for n ≥ 1. By setting f (η) = f 0 (η) + g(η) with f 0 (η) = η 3 /(η + L) 2 , the boundary conditions for f (η) are translated to corresponding conditions on g(η). In particular, f (η) = 1 as η → ∞ now becomes g (η) = 0 as η → ∞. Now, a solution of the form is sought by setting at η = η j = L 1+y j 1−y j for j = 0, 1, · · · , N, where {y j } N j=0 are the zeros of the Legendre polynomial of degree N + 1 (also known as the Gauss-Legendre nodes). Along with the conditions g(0) = γ and g (0) = 0, (56) leads to a system of N + 1 nonlinear algebraic equations in the expansion coefficients {a i } N i=0 in (55).

2.
Kajani et al. [59] truncate the infinite domain [0, ∞) to [0, L] for a "large enough" L and use what they call shifted Chebyshev polynomials to solve (5) and (6). The shifted Chebyshev polynomials are defined by in which {T n (y)} are the well-known Chebyshev polynomials. While the Chebyshev polynomials T n (y) are orthogonal on [−1, 1] with respect to the weight function σ(y) = 1/ 1 − y 2 , the shifted Chebyshev polynomials S n (η) are orthogonal on the interval [0, L] with respect to the weight function 1/ Lη − η 2 . The authors seek a solution of the form They use (58) to obtain in which the coefficients {b i } N+1 i=0 and {c i } N+2 i=0 may be easily expressed in terms of the coefficients {a i } N i=0 by evaluating the integrals needed to go from f (η) to f (η) and f (η). The solution to (5) and (6) is obtained by setting at η = η j = L 2 (y j + 1), for j = 1, 2, · · · , N, where the collocation points y j are the zeros of T N (y), and setting corresponding to the boundary condition as η → ∞. Note that the boundary conditions at η = 0 are automatically satisfied through the definitions of f (η) and f (η) as given in (59).

Other Methods
In this subsection, we will review a couple of methods that do not fall under the categories discussed previously in this paper.

1.
Using a truncated domain at η = η ∞ and imposing the asymptotic condition d 2 f /dη 2 = 0 as η → ∞, Asaithambi [60] solves (5) and (6) by formulating the problem in weak Galerkin form and with piecewise linear functions as basis functions and provides a finite-element solution. For this purpose, Asaithambi [60] rewrites (5) and (6) in the form in which the primes denote differentiation with respect to the variable ξ = η/η ∞ , g = f /η ∞ , and η = η ∞ is the truncated boundary which needs to be determined as part of the solution. Note that the semi-infinite domain [0, ∞) is now mapped to the finite domain [0, 1]. Subdividing the interval [0, 1] into N subintervals at ξ j = jh with h = 1/N and j = 0, 1, · · · , N, the basis functions φ j (ξ) defined by Asaithambi [60] constructs the solution in the form Thus, the method uses both finite-elements and finite-differences to obtain a nonlinear system of algebraic equations which is solved using Newton's method. However, this system has to be iteratively solved to adjust the value of η ∞ so that the asymptotic condition u (1) = 0 is satisfied. Asaithambi [60] handles this by using the secant method y setting and solving w(η ∞ ) = 0 for η ∞ . The superscript (k) here denotes the approximation to the different quantities during the k th iteration. The approximation η

2.
Recently, Hajimohammadi et al. [61] have developed a numerical learning approach to solving the Falkner-Skan problem (5) and (6). They cite several works from the literature that have used machine learning methods for the purpose of handling nonlinear models in the continuous domain, and demonstrate that a combination of numerical methods and machine learning methods could deal with nonlinearities much more easily and produce accurate solutions. In their work on the Falkner-Skan problem, they employ a combination of collocation using Rational Gegenbauer (RG) functions and the Least Squares (LS) Support Vector Machines (SM) and present results of their RG_LS_SVM method that compare well with results obtained previously by other researchers. The details of their method are too many to cover in this review and the reader is referred to the original publication [61].

Results
As mentioned previously in this paper, while existence and uniqueness results as well as results related to the properties of solutions have been explored extensively from theoretical perspectives, there is a significant amount of interest in exploring computational approaches for the solution of the Falkner-Skan problem. The researchers pursuing computational approaches have considered several special cases of (5) and (6) and demonstrated the effectiveness of their own approaches by comparing with those of others based on characteristics such as accuracy achieved, computational effort required, and ease of implementation. The flows described in Table 1 are most frequently studied in particular, and they correspond to the pairs of values (β 0 , β) in (7). The value of α = f (0) for each flow is also included in the table for reference. With β 0 = 1, most works [22,48,50] have reported results for both positive and negative values of β in the interval [−0.1988, 2]. These results are summarized in Table 2. Table 3 summarizes the results of previous works in which the authors have reported results for select values of β ∈ [10,40], with β 0 = 1.

Discussion
This article has provided a review of existing research on the analytical, semi-analytical, and computational methods for the solution of the nonlinear, third-order, two-point boundary-value problem known as the Falkner-Skan problem, as given by (5) and (6). The review includes brief descriptions of the approaches without going into the technical details. For analytical treatments, the extent of the analyses and the nature of the final results of the works considered have been presented. For numerical treatments, methods were discussed under three broad categories of shooting, finite-differences, collocation, and other (hybrid) methods. The results presented are representative of the most frequently reported results. The methods that have been described in some detail are intended to show the variety of approaches used to solve the Falkner-Skan problem. There are several other approaches to solving the Falkner-Skan problem that are closely related to those discussed in this review. Such approaches include differential transformation [62], variational iteration [63], iterative transformation [64] and other semi-analytical approaches [65] leading to series solutions.