A Review on a Class of Second Order Nonlinear Singular BVPs

: Several real-life problems are modeled by nonlinear singular differential equations. In this article, we study a class of nonlinear singular differential equations, explore its various aspects, and provide a detailed literature survey. Nonlinear singular differential equations are not easy to solve and their exact solution does not exist in most cases. Since the exact solution does not exist, it is natural to look for the existence of the analytical solution and numerical solution. In this survey, we focus on both aspects of nonlinear singular boundary value problems (SBVPs) and cover different analytical and numerical techniques which are developed to deal with a class of nonlinear singular differential equations − ( p ( x ) y (cid:48) ( x )) (cid:48) = q ( x ) f ( x , y , py (cid:48) ) for x ∈ ( 0, b ) , subject to suitable initial and boundary conditions. The monotone iterative technique has also been briefed as it gained a lot of attention during the last two decades and it has been merged with most of the other existing techniques. A list of SBVPs is also provided which will be of great help to researchers working in this area. [19] and extend the result of Shampine [21]. Finally, they demonstrated relative efﬁciencies of the various techniques on several real-life nonlinear SBVPs.


Introduction
The analysis of many important problems in science and engineering requires the solution of nonlinear singular differential equations. In the last few decades, several researchers explored the properties of solutions of nonlinear singular differential equations, and different methods like shooting method, nonlinear alternative, and upper-lower solution technique, etc. have been utilized. Difficulties arise when we deal with nonlinear singular boundary value problems (BVPs). The aim of this survey is to explore properties of nonlinear singular differential equations subject to certain initial and boundary conditions. First of all, we classify singularities and their occurrence in practical applications.
Consider the following general linear homogeneous differential equation of order 2 P(x)y (x) + Q(x)y (x) + R(x)y(x) = 0, where P(x), Q(x), and R(x) are continuous. A point on which we wish to study the properties of solutions of a differential equation (DE) may be classified depending on the behavior of P(x), Q(x), and R(x) around that point.
If P(x 0 ) = 0, x 0 is referred to as an ordinary point, and (1) is referred to as a regular differential equation at x = x 0 .
If P(x 0 ) = 0 but Q(x 0 ) = 0 or R(x 0 ) = 0, then x 0 is referred to as a singular point (SP) and (1) is referred to as singular differential equation at x = x 0 .
If P(x 0 ) = 0 and along with that (x − x 0 ) Q(x) P(x) and (x − x 0 ) 2 R(x) P(x) both are analytic at x 0 , then x 0 is referred to as a regular singular point; otherwise, it is an irregular singular point.
The SPs of a DE are usually few in number, so one can simply wish to ignore. However, these SPs often govern the principal features of the solution. Near an SP, a solution often becomes large or experiences rapid change or might be peculiar in some other manner. Thus, the behavior of a physical system modeled by a DE frequently is most interesting in the neighborhood of an SP. Often, geometric singularities in the modeling of the physical situation, such as corners or sharp edges, lead to SPs in the corresponding DE. Thus, precisely at these points, it is necessary to study the solution most carefully.
To characterize the behavior of the solutions of the DE (1) near x = x 0 , we require additional information about Q/P and R/P. Sometimes both linearly independent solutions are bounded as x → x 0 , or sometimes only one is bounded and the other is unbounded as x → x 0 , or sometimes both may be unbounded as x → x 0 .
Finally, if the DE (1) have solutions that become unbounded as x → x 0 , it is often important to determine how these solutions behave as x → x 0 . For example, does y(x) → ∞ in the same way as 1 (x−x 0 ) 1/2 or as 1 (x−x 0 ) or log (x − x 0 ) or in some other way? Therefore, to understand a physical process governed by such singular DEs, it is necessary to analyze their analytical properties especially in the neighborhood of an SP.
We arrange the paper in sections and subsections as follows. In Section 1.1, we present the role of boundary conditions (BCs). We also discuss the impact of these conditions on the singular differential equations. Furthermore, in Section 1.2, we show the origin of the different classes of SBVPs arising in different branches of sciences and engineering. We also provide a list of test examples that can be used as test cases. In Section 2.1, we describe numerical results and, in Section 2.2, we describe the analytical results. We have dedicated Section 2.3 to a monotone iterative technique, and its combination with other existing methods. Finally, in Section 3, we give future scope and, in Section 4, we conclude the survey.

The Boundary Conditions
Consider the following singular differential equation: − p(x)y (x) = q(x) f x, y, py , If f (x, y, py ) is continuous on [0, b] × D × D (D ⊆ R, D ⊆ R), p(x), q(x) is continuous on [0, b] and p(x) has no zeros on the interval [0, b], then the smoothness of solution of Equation (2) is at least C 2 [0, b]. However, if p(x) has a finite number of zeros in the interval [0, b], and f (x, y, py ) is continuous, then we can not easily make conclusions about the smoothness of the solution. If p(0) = 0 and p(x) is continuous, p(x) > 0 in (0, b], then, in the interval [0, b] for a small neighbourhood of 0, (2) behaves as a first order differential equation, and then it becomes difficult to predict about both solutions near zero. To maintain the regularity, various types of boundary conditions may be imposed on the DE (2). These boundary conditions can be decided by the behavior of p(x) and q(x) and f (x, y, py ).
In case f is continuous, p(0) = 0 and p > 0 ∈ (0, b] and b 0 ds p(s) = ∞ and conditions defined by (a * ), (b * ) are true, then, for any y ∈ D(L), we have y (0) = 0. In this situation, y ∈ C 2 [0, b]. In this case, we may consider y (0) = 0 or lim x→0 + p(x)y (x) = 0 at x = 0. We thus conclude the following: (d) If no singularity is present at the other end x = b, then the mixed type of boundary conditions may be imposed.
The physical and geometrical interpretations of a singular differential Equation (2) are well known. There are several real-life examples governed by such type of singular differential equations. In the next section, we discuss some physical phenomena arising in physiology, thermal explosion, astrophysics, etc. We describe some of the highly cited real life problems governed by DE (2).

Physiology
Oxygen Diffusion in a Spherical Cell: Lin [2] tried to predict the oxygen tension in a spherical cell by using an oxygen uptake kinetics of the Michaelis-Menten type. If V is the maximum reaction rate, P the oxygen tension and k m the Michaelis-Menten constant, then unsteady state oxygen diffusion in a spherical cell can be expressed mathematically as and the initial and boundary conditions are given by in which D is the diffusion coefficient of oxygen in the protoplasm, r 0 the radius of cell, h the permeability of membrane, r the radial co-ordinate, and t the time. For the convenience of computational purposes, Equation (3) and the initial and boundary conditions can be transformed into dimensionless form by introducing the following dimensionless variables and parameters:

Equations (3)-(6) are then transformed into
∂C ∂τ subject to McElwain [3] considered the steady state and reduced (7)-(10) to a singular ordinary differential equation subject to Anderson et al. [4] presented pointwise upper and lower bounds for the solutions of above SBVP.
Heat Conduction Model in a Human Head: The following boundary value problem arises in the study of the distribution of heat sources in a human head (see [5]) subject to Here, q(θ) is the heat production rate per unit volume, θ the absolute temperature, x the radial distance from the center, κ the average thermal conductivity inside the head, β a heat exchange coefficient, and θ a the ambient temperature. The form of q(θ) is the most important in this model. Flesh [6] assumed two possibilities on q as q = q 1 · θ 2 and q = q 2 · θ 3 , where q 1 and q 2 are constants. Gray [7] considered another approximation of q(θ) as where α and N are constants such that N is large subject to q > 0. He also observed that the law only holds over a limited range of temperature and would certainly not be applicable for temperatures in the hyperthermic region. Anderson et al. [8] generalized q(θ) by considering a nonlinearized form of q(θ) given by They also showed significance of the nonlinearity from a comparison with the linearized model. For smaller values of Nθ α , this q(θ) reduces to the form considered by Gray [7].

Electrohydrodynamics
Keller [9] expressed the equilibrium of a uniformly charged gas in a perfectly conducting container in terms of pressure p, the mass density ρ, and the charge density aρ (a is the ratio of electric charge density to mass density in the fluid) by the following equation: where is a non-negative increasing function of v = p p 0 dp ρ(p) . If we take p 0 = 1 and fluid is the ideal gas, then the differential Equation (16) becomes where Here, T is the constant temperature, R the gas constant, and m the average mass of the molecules in the gas. If the container is a sphere, a cylinder or a pair of parallel planes, it can be assumed that the solution u of the differential Equation (17) is a function of one variable only. This variable, which is denoted by r, is the distance from the center of the sphere, from the axis of the cylinder, or from the median plane in the three, two, or one-dimensional cases, respectively. If u = u(r) and n denotes the dimension, then the differential Equation (17) becomes The regularity of u at the center of sphere or axis of the cylinder requires that

Thermal Explosions
Thermal balance between heat generated and conducted away during a chemical reaction is given by where T is the temperature of gas, Q the heat of reaction, λ the thermal conductivity, W the reaction velocity, and ∇ 2 the Laplacian. Chamber [10] assumed that the reaction is monomolecular and its velocity follows the Arrhenius law: Here, k = 0 corresponds to the infinite plane-parallel plate vessel, k = 1 to the cylinder of length much greater than its radius, and k = 2 to the spherical vessel.
To express the temperature dependence of the rate constant, Nakamura et al. [11] proposed the following: where A, E 0 , and T 0 are parameters and n is an integer. Verma et al. [12] proposed and explored the following class of singular DE:

Polytropes as a Simple Model of a Star
The following derivation is due to Chandrashekhar [13]. Polytropes can be described by the equation of state in which the pressure P is given as a power-law where ρ is density and κ, γ are constants.

Epitaxial Growth
The Epitaxial Growth model is explored in detail by Carlos et al. in [14]. Let σ : Ω ⊂ R 2 × R + → R describe the height of the growing interface in the spatial point x ∈ Ω ⊂ R 2 at time t ∈ R + . σ satisfies the fourth order PDE where η(x, t) models the incoming mass entering the system through epitaxial deposition and λ measures the intensity of this flux.
The stationary counterpart of the partial differential Equation (33) subject to homogeneous Dirichlet boundary condition (35) and homogeneous Navier boundary condition (36) is defined as where η(x, t) ≡ G(x) is a stationary flux, and n is unit out drawn normal to ∂Ω. Let r = |x| and σ(x) = φ(|x|), and, due to symmetry, we arrive at where = d dr . If G(r) = 1, i.e., the new material is being deposited uniformly on the unit disc and, by using lim r→0 + rφ (r) = 0, w = rφ and integrating by parts from Equation (37), we have The differential Equation (40) is nonlinear, non-self-adjoint, and singular at the end r = 0.

Shallow Membrane Caps
Rachůnková et al. [15] explored the following SBVP which arises in the theory of shallow membrane caps (t 3 u (t)) + t 3 1 whereã,b, andc are given constants.

Catalyst Diffusion Reaction
Flockerzi et al. [16] propose the following system of SBVPs which model the class of catalyst diffusion reaction where L 1 , L 2 are linear expressions in y and β 1 , β 2 , K ij , i = 1, 2, j = 1, 2 are constant terms.

Remark 1.
It will be interesting to consider the problem arising in shallow membrane caps governed by Equations (41) and (42) because it contains the highly nonlinear term 1/u and 1/u 2 . One may also consider system on SBVPs of the type (43) and (44), as they are not explored much and investigations are still pending.

Test Examples
In this section, we present a list of test examples of the form (2). One may also refer a recent book by Agarwal et al. [17]. The problems for which the exact solution exist, we can modify to define general class of SBVPs where the nonlinear source term may depend on the derivative, i.e., f ≡ f (x, y, py ) and p(x) and q(x) may be some complicated functions: BCs: Exact Solution: .

Survey
Several numerical techniques, like finite difference method, collocation method, Galerkin method, shooting method, iterative method, finite element method, wavelet method, Lagrangian method, explicit direct method, Euler method, corrector-predictor method, extrapolation-extrapolation method, etc. and analytical techniques like topological degree method, monotone iterative technique, upper and lower solution method, etc. are available in the literature. To understand the development systematically, we divide this section into two subsections. In Section 2.1, we provide the important numerical results and in Section 2.2 we provide important analytical results. The monotone iterative technique has been described in Section 2.3 as in recent times a lot of attention has been paid to a monotone iterative technique.

Survey on Numerical Results
In the early forties of twentieth-century researchers suggested that the analytical solution of the singular boundary value problem (25)-(28) is possible only for k = 0, which they obtained. However, for k = 1 and k = 2, they employed numerical integration. Chamber [10] showed that analytical solution of the singular boundary value problem (25)-(28) for k = 1, 2 can be found and computed it for k = 1 in the terms of quadratures and for k = 2 in the terms of known tabulated function. Probably this was the first result which started motivating researchers to explore analytically further possibilities in the singular boundary value problems. Now, for convenience, we further divide the survey in two parts.

Numerical Results when Nonlinear Term is f (x, y)
Keller [9] considered the boundary value problem (18)- (19) and proved the existence of solutions. He used monotone iterative technique. In the late 1960s, Varga coauthored several papers on numerical techniques with Ciarlet, Natterer, Schultz, and others on regular as well as singular boundary value problems. They worked mainly in Sobolev Spaces. Ciarlet et al. [18] considered the following singular boundary value problem and suggested that, if 1 p ∈ L 1 [0, 1], then the singular boundary value problem (45)- (46) can be transformed into the following regular boundary value problem: by a simple change of variable Jamet [19] studied finite difference schemes rigorously for the following singular differential equation: where f (x) ∈ C(0, 1], f (x) → ∞ as x → 0, g(x), H(x) ∈ C[0, 1] and g(x) ≥ 0 with boundary conditions where σ is a number and 0 < σ < 1 or with the boundary condition He referred boundary value problem (50), (53) and (54) as one point boundary value problem since he required only u(x) to be bounded at x = 0 and observed that differential operator in (50) can be written as where p(x) = exp − 1 x f (t)dt . He discussed two finite difference schemes, a direct central difference analog of the Equation (50) and a scheme obtained by differencing the Equation (55). He proved that the differential Equation (50) has a unique solution and the solutions of both finite difference schemes converge to it. He also obtained a principal error function for the case f (x) = σ x , g(x) ≡ q (a constant) and illustrated a method of uniform error estimation. He showed the error to be O(h 1−σ ).
In the 1960s, Baily, Shampine, Waltman, Russell, and others studied regular nonlinear boundary value problems and later started analyzing singular nonlinear boundary value problems. In a classic work, Russell et al. [20] wrote that, if a physical law is expressed by and one is interested in planar, cylindrical or spherical geometries, he is led to the following singular differential equation: with k = 0, 1 or 2, respectively. They considered the differential equation (57) with the boundary conditions Using monotone iterative technique, they established existence-uniqueness results for the singular boundary value problem (57)-(58) for k = 1, when Lipschitz constant K > − j 0 b 2 (j 0 is a first positive zero of Bessel function of the first kind of order zero) and for k = 2, when K > − π b 2 . Furthermore, they examined three numerical techniques, namely collocation method, finite differences, and patch bases.

Remark 2.
These results modify the result of Jamet [19] and extend the result of Shampine [21]. Finally, they demonstrated relative efficiencies of the various techniques on several real-life nonlinear SBVPs.
In the 1970s and 1980s, M. M. Chawla established several results on the numerical solutions of singular boundary value problems. Consider the following singular boundary value problem: where A, B are finite constants, the source function f (x, y) is continuous, ∂ f /∂y exists, continuous In case b 0 ∈ (0, 1) based on three-point finite difference approximation, Chawla et al. [22] presented three methods for the singular boundary value problem (59) and (60) which are O(h 2 ) convergent. First, the method is on non-uniform mesh with one evaluation of f and the rest are on uniform mesh with three evaluations of f . Let us consider the following singular boundary value problem: where b 0 ≥ 0, a ∈ R. Excellent numerical treatments exploring various aspects of solutions of singular boundary value problems (59) and (61) with the boundary conditions (60) or (62) for different b 0 and ∂ f /∂y can be found in Chawla et al. [23][24][25][26][27][28], Katti et al. [29], Iyengar et al. [30], Sakai et al. [31,32], and Jain et al. [33] and the references therein. Now, we consider the following singular boundary value problem: Pandey [34][35][36] and Pandey et al. [37] generalized many results existing in the literature for a class of functions p(x) and established O(h 2 ) convergent methods for the singular differential Equation (63) with the boundary conditions (64) (α = 1, β = 0), y(0) = A or y (0) = 0, and The next result is due to Abu-Zaid et al. [38] in which they constructed an O(h 2 ) convergent difference scheme on a uniform mesh for the following singular differential equation: subject to the boundary conditions (py )(0 + ) = 0, y(1) = 0, where p(0) = 0 and p also satisfies several other hypotheses. Sen et al. [39] considered the following singular differential equation where one or more p k can be infinite at x = 0 with the boundary conditions and presented consistency, stability, and error estimates. The authors mainly extend the method of Gustafsson [40]. El-Gebeily et al. [41] considered the following singular differential equation: with the boundary conditions lim x→0 + p(x)y (x) = 0, y(1) = 0, where the functions w, p, q, and f : (0, 1) → R are Lebesgue measurable and include both limit point and limit circle cases. They gave a finite difference method under several assumptions on w, p, q and f . When w = p = x b 0 ≥ 0, the order of the convergence reduces to O(h 2 ) and verifies the methods given by other researchers. Guoqiang et al. [42] discussed a correction method and an extrapolation method for the following singular differential equation: subject to the boundary conditions u (0 + ) = 0, u(1) = A for b 0 ≥ 1. Pandey [43] established a second order spline method for the singular differential Equation (66) subject to the boundary conditions u (0) = 0 (b 0 ≥ 1) or u(0) = 0 (0 < b 0 < 1) and α 1 u(1) + β 1 u (1) = γ 1 . Now, consider the following singular differential equation: Ha et al. [44] discussed and analyzed numerical solutions of the above singular differential equation subject to the boundary conditions α 1 u(a) + β 1 u (a) = 0, α 2 u(b) + β 2 u (b) = 0, which they derived by the Green's function and shooting method. A numerical method consisting a simple iterative method coupled with a shooting method was presented by El-Gebeily et al. [45] for the differential equation Pandey et al. [46][47][48] consider the following singular differential equation with the boundary conditions and and extend the work of Chawla et al. [22], Chawla [23] for a class of functions p(x). The Adomian decomposition method (ADM) is discussed by Mittal et al. [49] for the singular differential Equation (67), where p(x) = x γ g(x), 0 ≤ γ < 1 and (∂ f /∂y) ≥ 0 with boundary condition (68). Pandey et al. [50] establish the O(h 4 ) convergent method for the singular differential Equation (63) with the boundary condition (69), which extends a result due to Chawla et al. [27] for a class of functions p(x).
The differential transform method is used by ASV Ravi Kanth et al. [51] for solving the following singular two-point boundary value problem, Bataineh et al. [52] use a modified homotopy analysis method (OHAM) and approximated solution of the following singular two-point boundary value problems: where p, q, r and g are continuous functions on (0, 1] and the parameters α 1 , α 2 , β are real constants. Ramos [53] presented series solutions of the Lane-Emden type equation subject to initial conditions based on either (i) formulation of Volterra integral equation or (ii) expansion of the dependent variable in the original ordinary differential equations.
It was shown that the series solutions were the same as those obtained by a direct application of ADM (Adomian [54]) to the original differential equation, He's homotopy perturbation technique (HPM)(He [55][56][57]), and Wazwaz's [58] two implementations of the ADM.

Remark 3.
ADM is an iterative method that depends on Adomian's polynomial. One may change the polynomial for better accuracy. To understand it, let us consider f (u) = u 2 . Now, we approximate u by ∑ u i and f (u) by ∑ A i , where A i 's are called Adomian's polynomial Adomian's polynomial omits the rest of the term. We can modify this approximation which may develop improved versions of ADM. Gebeily et al. [59] approximated the solutions of a class of nonlinear second-order singular boundary value problems (57)-(58) with a self-adjoint linear part. Taking advantage of certain boundary conditions, they obtained well-behaved solutions and, to avoid the singularity, they integrate and prove uniform convergence for the approximate solutions. For the various values of the deficiency index associated with the differential equation, they have shown the construction. A globally convergent iterative method is used for computation.
Caglar et al. [60] considered the differential Equation (57) subject to the boundary conditions defined by the Equation (69) and used a family of third-degree B-Spline to compute the solutions. They tested the method on several test problems.
Hasan et al. [61] developed a modified ADM to solve higher order nonlinear SBVPs defined as where N is a nonlinear differential operator of an order less than n, g(x) is given function and a 0 , a 1 , · · · , a n−1 , B are given constants. However, as test examples, the authors took only the case of n = 2, 3. Bataineh et al. [62] used homotopy analysis method (HAM) and approximated the solutions of singular initial value problems (IVPs) of the Emden-Fowler type given as

Remark 5.
Bataineh et al. [62] remark that HAM solutions contain an auxiliary parameter, which provides a convenient way of controlling the convergence region of the series solutions. They also have shown that the solutions obtained by the ADM and the HPM are only special cases of HAM solutions.
Ebaid [63] proposed an improved ADM to solve a class of linear and nonlinear singular boundary value problems including a Thomas Fermi equation subject to the Dirichlet boundary condition. Mohammadi et al. [64] used a Legendre operational matrix and computed the numerical solution of linear and nonlinear differential equations, for both initial and boundary singular value problems. It was based on the truncated Legendre wavelet expansion and collocation methods.
Shang et al. [65] presented HPM based on Green's function to solve (57)- (58). Its advantage to overcome the singularity behavior at the point x = 0 is discussed.
Secer et al. [66] applied the Sinc-Galerkin method to an approximate solution of second-order singular Dirichlet-type boundary value problems, which is based on approximating functions and their derivatives by using the Whittaker cardinal function. Remark 6. Secer et al. [66] have also shown an application of approximation of a functional involving nonlinear term which helps researchers to develop a new discrete scheme.
Bhrawy et al. [67] used a shifted Jacobi-Gauss collocation spectral method and solved the nonlinear Lane-Emden type equation where a, b 0 and b 1 are reals. The spatial approximation is based on shifted Jacobi polynomials P (α,β) T,n (x) with α, β ∈ (−1, ∞), T > 0, and n is the degree of the polynomial. The shifted Jacobi-Gauss points are used as collocation nodes.
Iqbal et al. [68] derived optimal OHAM for singular initial value Lane-Emden type problems and checked the effectiveness and performance. The method is simple and does not need discretization. OHAM provides a convenient way to control the convergence by optimally determining the auxiliary constants. It is remarked that this method converges rapidly at a lower order of approximations. Rismani et al. [69] proposed improved Legendre-spectral method to solve Lane-Emden type equations given as Rach et al. [70] considered the following system of Lane-Emden differential equations with given boundary conditions: The authors derived the modified recursion scheme for the components of the decomposition series solutions. The numerical results display that the ADM gives a reliable algorithm for analytic approximate solutions of these systems. The error analysis of the sequence of the analytic approximate solutions is performed by using the error remainder functions and the maximal error remainder parameters, which demonstrated an approximate exponential rate of convergence.

Remark 7.
Numerical results show that the method due to Babolian et al. [71] is superior (in terms of accuracy and speed of computing) to all related families of the spline method and the difficulties due to the singularity at x = 0 can be easily overcome here. Moreover, the numerical results are in excellent agreement with those obtained by the OHAM, variational iteration method (VIM), and the given families of spline methods.
Singh et al. [72][73][74][75][76][77] presented various methods based on ADM, VIM, and HAM for a class of singular and doubly singular BVPs. These results extend and generalize several results existing in the literature. Singh et al. [75] proposed two new modified recursive schemes for solving a class of doubly singular two-point boundary value problems defined as Here, a 1 , c 1 , a > 0, b ≥ 0, and c are finite constants. The condition p(0) = 0 ensures that the problem is singular and, in addition to this, if q is allowed to be discontinuous at 0, it is called doubly singular (Bobisud [78]). These schemes are based on ADM and integral operators. They avoid solving a sequence of nonlinear algebraic or transcendental equations for the undetermined coefficients with multiple roots, which is required to complete calculation of the solution by several earlier modified recursion schemes using the ADM. The approximate solution is computed in the form of series with easily calculable components. Singh et al. [76] considered (77) with q = p, subject to boundary conditions (68) or (64). They further assume that p(0) = 0. For y(0) = 0, they assume that They proposed a modification of the ADM. The technique depends on transforming the BVPs to Fredholm integral equations before establishing the recursive scheme for the solution components of a specific solution. The major advantage of the proposed method over the classical ADM or modified ADM is that it provides not only better numerical results, but also avoids unnecessary computation for determining the unknown parameters. Singh et al. [73] consider (57) subject to boundary condition (69) and developed an effective analytical technique based on the VIM in the presence of the convergence control parameter. To avoid solving a sequence of nonlinear algebraic or complicated integrals for the derivation of unknown constant, the boundary conditions are used before designing the recursive scheme for solution. The series solutions are found which converges rapidly to the exact solution. Convergence analysis and error bounds are discussed. Singh et al. [77] considered the problem of oxygen diffusion in a spherical cell with nonlinear oxygen uptake kinetics defined by (11) subject to C (0) = 0, C(1) + HC (1) = H. They transformed the BVP into an equivalent Fredholm integral equation and used the OHAM. Singh [79] solved (66) along with (69) and y (0) = 0, y (1) = c by a modified HPM. The advantage of the proposed schemes is that they do not require the computation of undetermined coefficients, whereas most of the previous recursive schemes do require the computation of undetermined coefficients. The author considered the integral form of the Lane-Emden equation and derived the recursive relation to manage the boundary conditions. The author also decomposed the domain into two subintervals for Neumann boundary conditions to deal with singularity. Singh [72] considered (77) subject to (78) and (79) and solved it by an iterative scheme. His method is based on Green's function and a modification of the HAM. Singh [74] considered (66) subject to (69) with (α 1 > 0 and α 1 = 0 (Neumann BC)) and the proposed algorithm is based on the HPM and integral form of the Lane-Emden equation. The proposed methods provide the direct recursive schemes for computing approximation to solutions. Remark 8. The advantage of these schemes by Singh et al. [72,73,[75][76][77] is that they do not require the computation of undetermined coefficients, whereas most of the previous recursive schemes do require the computation of undetermined coefficients.
Aydinlik et al. [80] use the Chebyshev finite difference method and analyzed the following class of SBVPs: where α ≥ 0, a 1 , a 2 , b 1 , b 2 , a and b are real constants and f (x, y) is a continuous real valued function and provided convergence and error analysis. Roul and others [81][82][83][84][85][86][87][88][89][90][91] presented numerical solutions for a class of SBVPs by using various methods e.g., HPM, decomposition HPM, optimal quartic B-spline collocation method (OQBSCM), high-order B-spline collocation method on a uniform mesh, the improved version of the standard HAM, the sixth-order numerical method based on sextic B-spline basis functions, the compact finite difference method (CFDM), the modified ADM and quintic B-spline collocation method, the B-spline collocation method, and the least square recursive approach (LSRA).
Roul et al. [90] considered (59) with p(x) = x α , subject to (68) and (69). They derived a new recursive scheme based on HPM and singularity was removed by transforming the original differential equation into an equivalent integral equation. Then, a recursive scheme without unknown constants was established by using HPM to approximate the solution. The convergence analysis is also given in the paper. Roul [92] convert (63) subject to (68) to an equivalent integral equation and employ the boundary condition to compute the undetermined coefficient. Roul [82] solved (63) subject to (68) and (69) by a modified HPM algorithm, called the domain decomposition HPM. The essence of the approach is to split the domain of the problem into some non-overlapping subdomains. In each subdomain, a method based on a combination of HPM and integral equation is implemented.
Roul et al. [88] considered the following SBVP: where α > 0, a > 0, and c are finite constants. They assume that g, ∂g ∂y both are continuous and ∂g ∂y ≥ 0, p(0) = 0 and p(x) > 0 in (0, 1]. They present a high-order numerical approach based on optimal quartic B-spline collocation method (OQBSCM) to solve the problem (82) and (83). Convergence analysis and global error estimates are also given. Thula et al. [91] used a high-order B-spline collocation method on a uniform mesh and solved (63)- (69). They claim that that their method has fourth-order convergence and is more accurate than finite difference methods, but their claim seems to not be fully justified (see [91], Table 9). Roul et al. [85] solved (63) subject to (68) and (69) by an iterative technique based on domain decomposition OHAM. This method produces approximate solution in the form of a series with easily computable components. Roul [83] presented a recursive scheme based on an improved version of the standard HAM to solve (66) along with (68) or (69). The author claims that like the FDM, FEM, or the spline method, HAM do not need linearization or discretization of variables and, like the perturbation methods, the present method does not require any small parameters to be present in the problem. HAM provides the solution in the form of a rapid convergent series with easily calculable components and provides a convenient way to control or adjust the convergence region of the series solution. Roul's HAM differs from classical HAM in the sense that it provides an optimal value for h, which greatly accelerates the convergence of the series solution.
Remark 9. The claim in the conclusion of the paper by Roul [83] that it does not require linearisation etc seems false and requires further investigation.
Roul et al. [87] considered the following second-order two-point singular boundary value problem governing electrohydrodynamic flow in a circular cylindrical conduit: where H is the Hartmann electric number and developed a sixth-order numerical method based on sextic B-spline basis functions.

Remark 10.
The technique of sixth order optimal B-spline collocation method can be further extended to solve general class of singular boundary value problem.
Madduri et al. [81] considered the following system of Lane-Emden equations: subject to the boundary conditions The authors present a fast-converging iterative scheme to approximate the solution by transforming the system into an equivalent system of Fredholm integral equations. The resulting system of integral equations is then efficiently solved by the optimized HAM.

Remark 11.
There are only a few papers on a system of SBVPs and most of the researchers focus on coupled systems with separated boundary conditions of the type y (0) = 0, y(1) = α, z (0) = 0, z(1) = β. One may use a coupled type of boundary conditions (y and z are not separated) and integral boundary condition to develop different numerical schemes.
Roul [84] consider (66) & (69) and proposed a numerical method based on a combination of modified ADM and quintic B-spline collocation method. The principal idea of this approach is to decompose the domain of the problem into two subdomains. In the first domain, the underlying singular boundary value problem is efficiently tackled by a modified ADM and, in the second domain, a collocation approach based on quintic B-spline basis function is designed for solving the resulting regular boundary value problem. Roul et al. [89] solved (66) and (69) (with β 1 = 0) by two B-spline collocation methods. The first method is based on a non-optimal quintic fourth-order convergent B-spline collocation approach, while the second method is based on the sixth order convergent optimal quintic B-spline collocation approach. Roul et al. [86] presented a technique, called the least square recursive approach (LSRA), to approximate the solution of the following Lane-Emden type initial value problems Niu et al. [93] considered (66)- (69) proposed an efficient method based on the simplified reproducing kernel method (SRKM) and the least squares method (LSM). They prove some new error estimations and also improve existing error estimations.

Remark 12.
We can also apply technique due to Niu et al. [93] to equations (84) and (85) and explore further.
Verma et al. [94] analyze Mickens' type non-standard finite difference schemes (NSFD) and establish their convergence and applied these schemes on Lane-Emden Equations (66)- (69). They compared numerical results with existing analytical solutions or with solutions computed with standard finite difference (FD) schemes. The result shows that the NSFD behaves qualitatively in the same way as the original equations. The NSFD approximate solution near a singular point efficiently where FD fails to do so (Buckmire [95]). The major advantage here is that, for only few number of spatial divisions, authors achieve high accuracy.

Remark 13.
NSFDs are going to be the most favored techniques in the near future as they consume very little time and provide highly accurate solutions for few spatial divisions.
Singh [96] considered the system of the Lane-Emden initial value problem defined by the following equations: They considered an equivalent Volterra integral form and applied OHAM. Convergence and error analysis of the proposed method are also discussed. Singh et al. [97] considered (86) for i = 1, 2, n = 1, 2 subject to boundary conditions w i (0) = 0, w i (1) = c i for i = 1, 2 and use the HAM with Green's function to extend their work on IVPs (86) and (87) to BVPs.
Marioara et al. [98] used a least squares differential quadrature method (LSDQM) on several applications of the standard Lane-Emden equation defined as where N > 0 and f (t, y) is a nonlinear function. The accuracy of the method is illustrated by a comparison with approximate solutions previously computed using methods like Boubaker polynomials expansion scheme (Boubaker et al. [99]), Hermite function collocation method (Parand et al. [100]), Rational Legendre pseudospectral approach (Parand et al. [101]), and Lagrangian interpolation method (Parand et al. [102]). Marinca et al. [103] solved a nonlinear singular system of polytrophic spheres of the Lane-Emden type arising in astrophysics by an optimal homotopy asymptotic method. Their procedure is combined with a rational expression which depends on the optimal HAM solution. The results obtained through their approximate analytical approach are compared with numerical results. Now, we discuss some results involving wavelets. Kaur et al. [104] considered the Lane-Emden equations of the type which occurs in the theory of stellar structure for the gravitational potential of a self gravitating, spherically symmetric polytropic fluid which models the thermal behavior of a spherical cloud of gas acting under the mutual attraction of its molecules and subject to the classical laws of thermodynamics. Here, α is related to geometry. They used quasi-linearization approximation and replacement of an unknown function through a truncated series of Haar wavelet series of the function.
Singh et al. [74,105] solved (66) along with (68) or (69) or the Neumann boundary condition by the Haar wavelet collocation method. To overcome the singular behavior at the origin, they first transform the Lane-Emden equation into an integral equation. The Haar wavelet collocation approach is used to convert the system into a system of nonlinear equations which are then solved by the Newton-Raphson method.
Verma et al. [106] solved (66) along with (69) by using a Haar wavelet coupled with a quasilinearization approach (HWQA). To show the accuracy of the HWQA, several examples are presented. Convergence of the proposed method is also established in this paper, which shows that the proposed method converges very fast. Verma et al. [12,107] proposed Haar, Legendre, Hermite, Chebyshev, Laguerre, and Gegenbauer wavelets collocation methods and solved (66) and (69). Verma et al. [108] presented the Haar wavelet collocation method and solved the coupled Lane-Emden equation: where t > 0, subject to initial values, boundary values, and four point boundary values: where n 1 , n 2 , v 1 , v 2 ∈ (0, 1) and k 1 ≥ 0, k 2 ≥ 0, ω 1 < 1, ω 2 < 1, δ 1 , δ 2 , γ 1 , γ 2 are real constants. Results are compared with exact solutions for IVPs and BVPs and for four point BVPs. Convergence of these methods is also established and found to be of second order. Rasanan et al. [109] solved (88) and (89) by using an artificial neural network framework subject to certain initial conditions. Verma et al. [110,111] considered the following non-self-adjoint, singular, nonlinear fourth order boundary value problem which arises in the theory of epitaxial growth where = d dr . The problem depends on the parameter λ and admits multiple solutions. Therefore, it is difficult to pick multiple solutions together by any discrete method like a finite difference method, finite element method, etc. Here, the authors propose a new technique based on HPM and VIM. The authors also study the convergence analysis of both iterative schemes in C 2 ([0, 1]). Earlier, Escudero et al. [14] proved existence and non-existence of SBVPs defined by (90).

Remark 14.
Verma et al. [110] proposed iterative techniques to capture multiple solutions of the non-self-adjoint differential equation, which is nonlinear and singular. Since the problem contains a non-self-adjoint operator, it is not easy to predict the solution and prove theoretical results. Hence a lot of investigations are pending for such nonlinear SBVPs.
Singh et al. [112] proposed VIM coupled with HPM to solve (66)- (69). An approximate solution is computed in the form of a series, which is very handy from a computational point of view. The accuracy of the proposed method is revealed by several real life examples. Singh et al. [113] presented a modified version of VIM and used Bessel and the modified Bessel functions to define a relaxation parameter to compute a solution of SBVP (66) and (69) and generalized the result of Ravikanth et al. [114].

Remark 15.
The major drawback of the modified scheme proposed by Singh et al. [113] is that they could not compute the optimal value of the parameter.
Verma et al. [115] proposed fuzzy transformed finite difference methods (FTM) to solve and allowed the source term to have jump discontinuity. FTM solutions are in good agreement with exact solutions. In case of jump discontinuity in the source term, the authors observe that, with an increase in jump, the FTM method gives better results than the FD method.
Recently, Verma et al. [116] studied the Taylor series solution to find the exact as well as approximate solution of SBVP and also showed that this technique is more efficient than VIM as well as ADM.
Recently, Julee et al. [117] proposed an efficient numerical technique for the numerical solution of Lane-Emden-Flower type boundary value problems. They applied a Bernstein collocation method on an equivalent integral equation.

Remark 16.
It will be interesting to see the effect of fuzzy transforms on other existing techniques and their effect on nonlinear SBVPs [17].
Matlab Package: Ewa Weinmueller [118] developed a MATLAB package known as SBVP 1.0 (SingularBVP) to solve the first order BVPs of the form The function f on the right-hand side of the above equation may contain singularities of the first kind, f (t, y) = M(t) (t−a) + g(t, y), where M(t) is a matrix that depends continuously on t and g is a smooth vector-valued function. They used the method of collocation. The default choice is collocation at equidistant points, Gaussian points are optional. The code is suitable for a wide range of tolerances. For efficiency, the order of the method is chosen with respect to the prescribed solution accuracy.
Mathematica Package: Cabada et al. [119,120] developed the Mathematica package which provides a tool for calculating the explicit expression of the Green's function related to an nth order linear ordinary differential equation, with constant coefficients, coupled with two-point linear boundary conditions. Remark 17. Many years have passed since people started solving SBVPs and our computational facility has also seen tremendous improvement in terms of space, speed, and efficiency, but we still don't have a suitable package to solve a larger class of SBVPs. Thus, it will be interesting to see if we can develop Mathematica/Matlab packages that can solve most of the SBVPs.

Numerical Results when the Nonlinear Term is f (x, y, py )
Jain [121] has obtained one parameter family of a self starting Runge-Kutta-Nyström method for the solution of a general second order singular initial value problem with the spherical symmetry of the form where b 0 = 0, 1 and 2 subject to initial conditions u(0) = A, u (0) = 0. The methods are exact for u(x) = x −1 , 1, x, x 2 , x 3 and x 4 . Qu et al. [122] considered the following singular nonlinear two-point boundary value problem: y (0) = 0, y(1) = y r , and presented an iterative algorithm. Error estimates for uniform partitions are investigated, and it has been shown that for sufficiently smooth solutions the method produces O(h 4 ) approximations.
Roul et al. [123] takes the following class of BVPs: where µ > 0, η ≥ 0 and ρ is a finite constant. The authors claim that they have constructed and attained convergence of two B-spline collocation methods for the above class of nonlinear derivative dependent singular boundary value problems. Singh et al. [124] considered the following class of singular BVPs: (p(t)y (t)) = q(t) f (t, y, py ), where a > 0, b, and c are any finite real constants, p(0) = 0 and q are allowed to be discontinuous at x = 0. Furthermore, they assume that (i) p ∈ C[0, 1] ∩ C 1 (0, 1], p > 0 in (0, 1]; (ii) q > 0 ∈ (0, 1], q ∈ L 1 (0, 1] and q is not identically zero; (iii) x 0 q(s)dsdx < ∞; (iv) f (t, y, py ) is continuous, not identically zero and Lipschitz. They presented an algorithm based on the ADM and Green's function and proved the convergence. The accuracy is measured using the maximum absolute error or absolute residual error.
Singh et al. [125] consider (96) with Dirichlet boundary condition at x = 0, such that and extended the work of Singh et al. [124]. Roul [126] claims to derive that a fast-converging recursive scheme is presented to approximate the solution of a class of derivative dependent doubly singular boundary value problems (94) subject to (97) or (98). The original problem is reformulated as an equivalent integral equation and then was solved by an improved HAM with a parameter that accelerates the convergence. Roul et al. [127] developed a fourth order convergent optimal quintic B-spline collocation (OQBSC) method along with its convergence for the set of Equations (94) and (95). Roul et al. [128] developed a fourth order convergent compact finite difference method (CFDM) for solving a general class of two-point nonlinear singular boundary value problems with Neumann and Robin boundary conditions (94) and (95).

Remark 18.
There are some serious issues in assumption for the existence and uniqueness of the solution for the derivative dependent problem in most of the above discussed numerical methods. The authors do not mention the assumption on the Nagumo Condition that is required to control the derivative. In addition, assumptions by Roul et al. [123] on Equation (94) are not appropriate as in a derivative dependent case mostly authors always consider the nonlinear term as f (x, y, py ) not as f (x, y, y ). As the former one helps us in controlling the derivative, it is very difficult to have control over the derivative.

Survey on Analytical Results
This section is further divided into three subsections. First, we give a brief account of the techniques used so far. Then, we write analytical results when the nonlinear term is derivative independent. Later, we provide a survey when the nonlinear term is derivative dependent.

Some Existing Analytical Techniques
Several methods are available to study the analytical properties of singular boundary value problems. Broadly, they can be divided into the following three categories. This division is due to Zhang [129].

Shooting Method
The Shooting method had been used successfully to study some special singular problems like Negative Exponent Emden-Fowler boundary value problems (Taliaferro [130]). In fact, this method works well provided f (x, y) is decreasing in y. However, for other cases, it often seems useless.

Nonlinear Alternative
Leray et al. [131] in their celebrated paper introduced some "Nonlinear Alternative" theorems for compact maps. These theorems have enhanced greatly the theory of ordinary differential equations. There are two major approaches to modern nonlinear alternative theory, (i) Topological Degree Method: This is based on degree theory (Lloyd [132]). (ii) Topological Transversality: This is based on essential maps (Granas [133]).
Topological Methods, though it has many advantages in treating non-singular problems, still has some difficulties when treating singular problems.

Upper and Lower Solution Methods
Upper and lower solution techniques were developed originally for non-singular problems. Since there are many differences between singular and non-singular problems, using it as a general method to the singular problems needs basic investigations Recently, there has been a lot of activity related to the theory of upper and lower solutions. It has been successfully coupled with other existing techniques, e.g., topological degree theory (Duhox [134] [138]) and fixed point theory (Marcelli et al. [139], Calamai et al. [140], Biagi et al. [141]) to study the analytical properties of singular boundary value problems. Zhang [129] shows that the upper and lower solution technique is the most promising technique as far as singular boundary value problems are concerned. Zhang's claim is further justified by Pandey and Verma in their work [142][143][144][145][146][147]. Cabada et al. [148] discussed contemporary and classical results related with the existence of solutions of second-order problems coupled with nonlinear boundary value conditions.

Analytical Results when the Nonlinear Term is f (x, y)
Now, we discuss analytical results for the following class of singular differential equation: We have already discussed some analytical results for the singular boundary value problems due to Chamber [10], Keller [9] and Russell et al. [20].
Dunninger et al. [149] considered the following class of singular differential equation: with u(1) = 0, where p(0) = 0, p(x) > 0 in (0, 1] and p ∈ C[0, 1] ∩ C 1 (0, 1]. Boundary condition at the singular end is decided by the behavior of p(x) in the neighborhood of zero. In this work, they first derived priori bounds in an appropriate weighted Sobolev norm, which in turn imply L ∞ a priori bounds. Then, they applied the fixed point mapping of a cone to prove existence of positive solutions to the differential Equation (100). The next important result is due to Chawla et al. [150] in which they considered the singular differential Equation (61) with the boundary conditions for b 0 ≥ 1. They also assumed y ∈ C 2 (0, 1], f (x, y) to be continuous on S = {0 ≤ x ≤ 1, −∞ < y < ∞}, and ∂ f /∂y exists and continuous on S. They observed that, in most of the existence results available until then in the literature, it had been invariably assumed that −∞ < (∂ f /∂y) ≤ 0 (independent of b 0 ). Only for two cases b 0 = 1 and b 0 = 2, Russell et al. [20] proved the existence of unique solution for b 0 = 1 if Lipschitz constant K < j 2 0 (j 0 is the smallest positive zero of J 0 (x)) and for b 0 = 2 if K < π 2 . They therefore raised a question that, for the general nonlinear singular two-point boundary value problem (61)  Pandey [151][152][153] studied the singular differential Equation (99) extensively, where q(x) = p(x) and p(x) satisfy the following conditions: Pandey [151,152] considered both types of boundary conditions y(0) = A or y (0) = 0 and y(b) = B for b 0 dt p(t) < ∞. Ciarlet et al. [18] suggested that, in this case ( b 0 dt p(t) < ∞), the singular boundary value problem can be transformed into a regular boundary value problem but Pandey [151,152] proved that the direct consideration of the singular nonlinear boundary value problem gives better result. When where v 0 and u 0 are lower and upper solutions satisfying certain inequalities. These results improve the result of Ciarlett et al. [18] and extend the results of Russell et al. [20] and Chawla et al. [150] for a class of functions p(x). To establish these results for a class of functions p(x), series expansion of an arbitrary function in terms of orthonormal eigenfunctions is required. Thus, first Pandey [151][152][153] established three fundamental theorems, which are similar to Theorem 2.7 (i), (ii) and Theorem 2.17 of Titchmarsh [154] then used it to establish the existence and uniqueness results.
Remark 20. Pandey [151,152] proved that the direct consideration (without reduction as suggested by Ciarlet et al. [18]) of the singular nonlinear boundary value problem gives a better result, in the sense that it enlarges the class of functions.
A rigorous study of the analytical properties of the singular boundary value problems was done by O'Regan [155]. O'Regan [156] considered the following class of singular differential equation: with appropriate boundary conditions which can model both Thomas-Fermi and Emden-Fowler equations. In this work, the author studied the existence of solutions to the differential Equation (102) with the boundary conditions and Assuming further conditions on p and ψ, and using nonlinear alternatives of Leray-Schauder due to Granas et al. [157], he obtained various existence results. One important observation in this result is the sign restriction imposed on f (x, y) which says that there is constant M 0 ≥ 0 such that Existence and uniqueness of solutions of the singular differential Equation (99) on 0 < x < 1 with the boundary conditions lim x→0 + y(x) cos α + p(x)y (x) sin α = 0 and y(1) cos β + p(1)y (1) sin β = 0 was established by El-Gebeily et al. [158]. They first derived some results on the spectrum of an appropriate singular differential operator and then used it to show that, for certain boundary conditions, the problem has a unique solution whenever the range of (∂ f /∂y) has an empty intersection with the closure of the spectrum of the singular differential operator, where the nonlinearity f satisfies a "nonresonance" type condition. This work generalizes the work of Chawla et al. [150].

Remark 23.
Zhang [129] shows that the upper and lower solution technique is the most promising technique to tackle the singular boundary value problems.
The technique of Mahwin et al. [161,162] was used by O'Regan [160] to prove existence results for the following non-resonant singular boundary value problem of limit circle type where p ∈ C[0, 1] ∩ C 1 (0, 1), p > 0 on (0, 1), q is measurable with q > 0, a.e., on (0, 1] and Agarwal et al. [163] presented some existence results which complement and extend several existing results in the literature for the following singular boundary value problems: Conditions assumed on the nonlinear term f and q are not mentioned here due to brevity. They also considered the following singular boundary value problem: where p ∈ C[0, 1] ∩ C 1 (0, 1), p > 0 on (0, 1) and 1 0 ds p(s) < ∞, and suggested that Liouville transformation can be used to reduce the singular boundary value problem (112) and (113) to boundary value problem (110) and hence these ideas can be used to establish existence results for the singular boundary value problem (112) and (113).
Considering the singular differential Equation (110) with q = 1, where f (t, y) is not a Carathéodory function due to the singular behavior of its variables (t, y) and using the technique of Bobisud et al. [164], O'Regan [165] generalized and complimented the results of Habets et al. [166]. Nonlinear boundary conditions y(0) = 0 = θ(y (1)) + y(1) were also considered. Agarwal et al. [167] presented a new existence result for the singular differential Equation (99) (p = q = 1) with Dirichlet or mixed boundary data and allowed the nonlinearity to change sign with singularities at y = 0, t = 0, and t = 1.
In a different approach via topological degree theory coupled with upper and lower solutions, Duhoux [134] established the existence results for the following nonlinear singular Sturm-Liouville problem −(r(t)u (t)) + p(t)u(t) = f (t, u) over the compact interval [a, b], where f is a Carathéodory function. Ravi P. Agarwal presented several existence results that extend and complement many results in the literature for the boundary value problem (110) with a nonlinear term f of sign-changing type. For y(1) = 0, Agarwal et al. [136] used a theory of upper-lower solutions coupled with a Schauder's fixed point theorem, and for the case θ(y (1)) + y(1) = 0 they used a Nonlinear Alternative of Leray-Schauder type with a growth restriction on f . In case nonlinear term f may not be a Carathéodory function because of the singular behavior of the variable y, f may be singular at y = 0 was considered by Agarwal et al. [168,169]. The latter result relied on the well known Rayleigh-Ritz inequality.
Agarwal et al. [170] considered following class of singular BVP: where for some , C > 0 and γ ∈ (0, 1), so that it may be singular at y = 0 (C may depend on ), and used variational set up to obtain positive solutions.
Using an approximation method together with the theory of upper and lower solutions, Lü et al. [171] studied the existence of solutions to the singular differential equation −u (t) = g(t, u) + h(t, u), t ∈ (0, 1) subject to the boundary conditions u(0) = 0 = u(1), where a nonlinear term may be singular at t = 0, 1; u = 0 and moreover may change sign. El-Gebeily et al. [172] considered the singular differential Equation (99), a.e., on I = (0, 1) and established existence of solutions provided the problem has a lower solution and an upper solution. They use a quasilinearization method to obtain two monotonic sequences of approximate solutions converging quadratically to a solution of the differential Equation (99).
Ford et al. [178] studied (66) subject to (68)-(69) and proved existence and uniqueness of solutions. The proof is constructive in nature and could be used for the numerical generation of the solution. f (x, y) is not singular in x and singularity in y is easily avoided. Solutions are found in finite regions where ∂ f ∂y ≥ 0, using an integral equation whose Green's function contains an adjustable parameter that secures convergence of the Picard iterative sequence. Several well-known test examples are tested by their method.
Benmezai et al. [179] proved existence of positive solutions by using fixed point theory: where a ∈ C 1 ((0, 1), (0, ∞)), 1/a is integrable on any compact subset of (0, 1], b ∈ C((0, 1), [0, +∞)) does not vanish identically and is integrable on any compact subset of [0, 1), and f : [0, 1] × R + → R + is continuous with f (t, u) > 0 for all (t, u) ∈ [0, 1](0, ∞). Some test examples are also discussed. Fewster-Young et al. [180] generated new a priori bounds on the solutions to the second-order coupled SBVPs formulation and application of differential inequalities, −αy(0) where y ∈ R n and f is a continuous function from [0, T] × R n to R n . Based on shooting method given positive integer N, Baxley et al. [181] derived conditions on the nonlinear function f which guarantee that the boundary value problem y = f (t, y), 0 < t < 1, y(0) = 0, y (r) = 0 has N positive solutions. The nonlinear function is allowed to be singular at y = 0 and t = 0 but is required to satisfy an integrability condition reminiscent of a condition first used by Taliaferro [130] and growth conditions are similar to those assumed in the nonsingular case by Henderson et al. [182].
Escudero et al. [14] and Verma et al. [183] considered (90) a non-self-adjoint, singular, nonlinear fourth-order boundary value problem, which arises in the theory of epitaxial growth and proved existence and non-existence.
Afgan [184] considered a class of second order initial value problems y (t) + p(t)y (t) + q(t, y) = 0, y(0) = a, y (0) = 0, and extended the family of second order initial value problems for which the existence and uniqueness of a solution is known. To handle the singularity, they introduce the more general Carathéodory type of restrictions. He et al. [185] considered the following class of system of Lane-Emden equations and applied Taylor series to compute the exact solution defined as u(x) = [186] considered where p and q : (a, b) → R are locally integrable functions with singularities at the points a and b, u(a+) is the right limit of the function u at the point a, and u(b−) and u (b−) are the left limits of the functions u and u at the point b.
Kiguradze [187] established new tests for the solvability and unique solvability of two-point boundary value problems for ordinary second-order differential equations with non-integrable singularities in the time variable. The author also described a set of functions f : is satisfied for arbitrary x ∈ R and l > 0, but the boundary value problem u = f (t, u), u(a+) = 0, u(b+) = 0, has a unique solution.

Analytical Results when the Nonlinear Term is f (x, y, py )
Bernstein [188] considered the regular differential equation and established a sufficient condition for unique solvability subject to Dirichlet boundary conditions. Granas et al. [189,190] generalized the result of Bernstein [188].
Dunninger et al. [1] considered the following class of singular boundary value problems: where p, q ∈ C[0, 1] ∩ C 1 (0, 1], p(0) = 0 and p, q > 0 on (0, 1]. Boundary condition at t = 0 is decided by the behavior of p(x) in the neighborhood of t = 0. Nonlinear term f (t, u, pu ) is continuous on Sign Restriction: u f (t, u, 0) > 0 for |u| > M 0 and Nagumo Condition: where M 0 is a constant. Their main purpose was to extend the main result of Granas et al. [189,190] for the singular problems (123).

Remark 25. Nagumo condition is required to keep the derivative bounded.
A more singular but disjoint from Dunninger et al. [1] was considered by Bobisud et al. [191] given by with the boundary conditions y(0) = a and y(1) = b, where φ and ψ are continuous, ψ > 0 on exists and 1 0 tdt ψ(t) < ∞. Furthermore, they assumed that there is a constant M ≥ 0 such that y φ(t, y) > 0 for |y| > M. This result again extends the work of Granas et al. [192,193]. These were the initial results based on topological transversality, but there are many more in the sequel. Bobisud et al. [164] established existence of positive solutions of the following singular differential equation: In the two classes of problems, they considered f to be positive at least for small values of y and y ; in the first class, f ≥ 0 holds for all y and in the second f changes sign. They obtained a priori estimates on the solutions to related 2-parameter families of problems and then used these bounds with the topological transversality theorem and the Ascoli-Arzela theorem to get existence for the original problem. Based on a topological transversality theorem, O'Regan [194] improved the result of Bobisud et al. [164]. In another work, O'Regan [195] further analyzed the unexplored aspects of the singular differential Equation (125) and improved the results of Bobisud et al. [164] and O'Regan [194] (and the results given in the references), the tool was same, i.e., topological transversality theorem.
Bobisud [78] considered the following class of singular boundary value problems: (p(t)y (t)) = q(t) f (t, y, py ) on (0, 1), lim t→0 p(t)y (t) = 0, and established the existence of solutions under a variety of conditions. Since p(0) = 0 and q(x) is allowed to have integrable singularity, he referred to differential Equation (126) as 'doubly singular'. He also examined a singular differential Equation (126) subject to the Dirichlet boundary conditions y(0) = a and y(1) = b. He pointed out that the condition f (x, y, 0) having the same sign of y for large y is unnecessarily restrictive as the simple example y = 2 shows. For the Dirichlet boundary conditions, he gave an alternative Lemma free from sign restriction. This result again relied on topological transversality.
He used the upper and lower solution method to prove the results. This result improves the result of Dunninger et al. [1] and partially improves the result due to Bobisud [78] (Theorem 1). Zhang [177] without any real change allowed q(t) to be discontinuous at t = 0 and extended these results for the doubly singular problems.
The existence of non-negative solutions of the second-order singular differential Equation (128) subject to the homogeneous Dirichlet, Neumann, Robin or periodic boundary conditions was studied by Bobisud et al. [135]. The differential equation may be singular at t = 0 or t = 1. The method of lower and upper solutions is developed by using topological transversality for Dirichlet and Robin boundary conditions and coincidence degree for Neumann and periodic conditions. The Bernstein-Nagumo type growth restrictions are imposed on f . Since p and q are assumed to be positive, the conditions f (t, 0, 0) ≤ 0, f (t, C, 0) ≥ 0, where C is a positive constant imply that 0 is a lower solution and C is an upper solution, thus there exists a solution between 0 and C, which is non-negative. Their final section deals with the singular Dirichlet problem for the equation where µ 0 > 0 is the smallest eigenvalue of the associated homogeneous problem. Theorems establishing the existence of non-negative solutions are proved under various growth restrictions. The nonlinear alternative is used here. It is interesting that in this section f is required to grow sufficiently fast in the second variable, whereas most other results are proved for f growing no more than linearly in the second variable. Now, consider the following singular differential equation: (p(t)y (t)) = f (t, y, py ), a.e., on [0, 1].
O'Regan [197] studied the differential Equation (136) subject to Sturm-Liouville, Neumann, Periodic, or Bohr type boundary conditions, where f is a L 1 − Carathéodory function. In this paper, he utilized properties of the Green's function, Sturm-Liouville eigenvalue problem, and nonlinear alternative of Leray-Schauder type to establish several existence results. Another result for the singular differential equation satisfying either the boundary conditions lim t→0 + p(t)y (t) = 0 = y(1) or the Dirichlet condition y(0) = 0 = y(1) is due to Agarwal et al. [198]. Agarwal et al. [198] proved existence theorems using a fixed point theorem due to O'Regan [155]. These results generalize and complement several results given in the references of the paper.

Remark 27.
Agarwal et al. [198] called the differential Equation (137) singular if f is singular at y = 0; otherwise, non-singular (in either case f may be singular at t = 0 and/or t = 1).
Let us consider the following singular boundary value problem: where the nonlinear term f may be singular at (i) y = 0 and y = 0, (ii) y = 0 but not y = 0 and (iii) y = 0 but not y = 0. Agarwal et al. [198] examined case (ii) in its full generality, but case (i) and case (iii) received very little attention in the literature so Agarwal et al. [199] tried to fill this gap in the literature. O'Regan et al. [200] considered singular differential Equation (137) subject to the boundary conditions lim t→0 + p(t)y (t) = 0 = y(1), where f is allowed to change sign and 1 p may not be necessarily in L 1 [0, b]. In addition, f may not be a Carathéodory function because of the singular behavior of the y variable. Under suitable growth restrictions on f and the existence of certain upper and lower solution type functions, they proved the existence of a solution lying between lower and upper solutions. The proofs were based on an approximating method of a sequence of non-singular problems, the Schauder fixed-point theorem, and the Arzelà-Ascoli theorem. They also extended the idea for Dirichlet boundary value problems. Using the method of upper and lower solutions combined with properties of the consequent mapping, Ntouyas et al. [201] studied singular differential Equation (137) and pointed that Agarwal et al. [198] exploited the zero set of the nonlinearity f and obtained a priori bounds for the solutions (and their derivatives) without (in some way) growth restriction on the third variable of f .
Agarwal et al. [202] considered the singular differential Equation (125) with boundary conditions y(0) = 0, y(1) = 0, and singular differential equation (p(t)y (t)) + f (t, y, py ) = 0 with boundary conditions lim t→0 + p(t)y (t) = 0, y(1) = 0, where f (t, y, py ) is allowed to change sign and f may be singular at t = 0. They proved the solvability by a combination of the regularity and sequential techniques and the method of lower and upper functions. Then, using the Schauder fixed-point theorem and Arzela-Ascoli theorem, they established the complete solvability of the boundary value problems. This paper generalizes assumptions of some of the previous works given in the references therein. Agarwal et al. [203] considered the singular differential Equation (125) with the boundary conditions y(0) = 0 and y(1) = 0 or y (1) = 0. The nonlinear term f may be singular in the independent variable y = 0 for the boundary condition y(1) = 0 and in the independent variable y = 0 and y = 0 for the boundary condition y (1) = 0. They presented a new approach based on Schauder's fixed point theorem and the results extend and complement those in the literature. Agarwal et al. [204] considered the existence of a solution of the singular differential Equation (125) with boundary conditions y(0) = y(1) = 0, by the upper and lower solution method. The case when f is not dependent on y has been studied extensively in the literature, while the general case considered here has been given less attention. In this paper, not only an existence theorem is established where f : (0, 1) × (0, +∞) × R → R is continuous and may be singular at y = 0, but a method is presented for the construction of upper and lower solutions by using the truncation technique. An example shows that the existence theorem is new even if f is independent of the third variable. Xie et al. [205] considered the singular differential Equation (136) with boundary conditions lim t→0 + p(t)u (t) = 0 = u(1).
They assumed that the nonlinearity f may change sign and can be singular at u = 0 and t = 1 and presented sufficient conditions to ensure the existence of positive solutions. The proof is based on the method of approximation and upper-lower solutions. Yan et al. [206] considered the singular second-order differential Equation (138) with boundary conditions y(0) = 0, y (1) = 0, where f may be singular at y = 0 and y = 0. Using fixed point index theory, they proved existence of multiple positive solutions.
Yan et al. [207] studied the following class of SBVPs, on infinite interval by using fixed point index, where f may be singular at x = 0 and px = 0, and f (t, x, z) ∈ C(R + × R + 0 × R + 0 , R + ), here R + = [0, +∞), R + 0 = (0, +∞). O'Regan et al. [138] established existence theorems using quasilinearization in the presence of upper and lower solutions for general nonlinear second order singular ordinary differential equations is defined for almost all t ∈ I = (a, b) with −∞ ≤ a < b ≤ ∞ and F is a nonlinear operator. They assume 1/p, q, w ∈ L loc (I), w > 0 almost everywhere in I. Rachunkova et al. [208] discussed the properties of differential equation where a ∈ R \ 0, and f satisfies L p -Carthéodory condition on (0, T] × R 2 for some p > 1. Asymptotic behavior as t → 0 + is given as well as necessary-sufficient condition, which ensures that solution is in C 1 [0, T]. Pandey et al. [145][146][147]209] generalized their earlier result (Pandey et al. [142][143][144]) to the following class of SBVPs with derivative dependent source function. These results also compliment and generalize several other results existing in literature: Pandey et al. [145] considered (142) subject to the boundary condition They assume the following conditions on p(x), q(x) and f (x, y, py ).
Furthermore, they assume b 0 1 p(x) x 0 q(s) ds and the homogeneous boundary value problem has only a trivial solution. Here, they do not assume that x = 0 is a regular singular point as assumed by Pandey et al. [152,153]. In case q(x) = p(x) = x α , the condition (144) says that 0 ≤ α < 3. Pandey et al. [146] prove the existence of the solutions without uniqueness and complemented the results of Pandey et al. [145] when x = 0 is a regular singular point. Though uniqueness could not be achieved, but existence is proved for all α ≥ 0. The results by Pandey et al. [145,146] improve results due to Dunninger et al. [1], Bobisud [78] and Zhang [177] as we do not require sign restrictions on the nonlinear term. Furthermore, these works also generalize the work of Ford et al. [178].
Pandey et al. [147] considered (142) subject to the following boundary consitions and generalize the result by Pandey [151]. In case b 0 dt p(t) < ∞, the researchers (Agarwal et al. [163], Ciarlet et al. [18], Zhang [176,177]) have suggested to reduce the singular problem to regular one by a change of variable. However, Pandey [151] suggested that a direct consideration of singular problem provides better results. In this work, the authors consider the singular boundary value problem directly. In addition, here they do not assume that the point x = 0 is of regular singular type as assumed in O'Regan [156] (Theorem 2.3).
In addition, in several results sign restrictions (i) y f (x, y, 0) > 0, |y| > M 0 , M 0 is a constant (Dunniger et al. [1], Bobisud [78], O'Regan [156] [176,177]) have been imposed. However, such sign restrictions are quite restrictive as the simple differential equation y = 2 does not satisfy the sign restriction (i) (Bobisud [78]). Pandey [147] does not require any sign restriction, so their results improve and generalize several results as mentioned. The authors assume that p(x), q(x), and f (x, y, py ) satisfy the following conditions: Furthermore, they assume that k(x) and µ(x) satisfy the following properties.
Finally, they assume that the homogeneous boundary value problem has only trivial solution. Pandey et al. [209] extend the result due to Pandey et al. [142] for derivative dependent source function, i.e., f = f (x, y, py ). Here, again sign restriction is not required. Now, we move on to monotone iterative techniques.

Survey on Monotone Iterative Techniques
Basic theory and common features of the monotone iterative technique for certain problems of ordinary and partial differential equations can be found in a book due to Ladde, Lakshmikantham, and Vatsla [210]. Historical development discussed below regarding upper and lower solution techniques can be found in a recent book by Coster and Habtes [211]. Here, we discuss briefly about upper and lower solutions related to a monotone iterative technique. If one wishes to work in this area we recommend the book by Coster and Habtes [211].
Actually, monotone iterative theory related to upper and lower solution goes back at least to E. Picard [212] who considered the Dirichlet problem and generated a convergent sequence of approximation (α n ) n from the constructive scheme −α n (t) + f (t, α n−1 ) = 0, α n (a) = α n (b) = 0 under the assumption that f is continuous function and decreasing in u. Under some more assumptions, he proved the existence of a positive function α 0 such that −α 0 + f (t, α 0 ) > 0 on (a, b) and α 0 (a) = α 0 (b) = 0. Such a function is called lower solution and this method is referred to as a monotone iterative technique. Dragoni [213] considered the Dirichlet boundary value problem assumed the existence of functions α and and obtained existence of a solution u together with its localization α ≤ u ≤ β. The regularity assumptions were that f is continuous and bounded on These results are probably the first ones that recognize explicitly the central role of the functions, α and β, which are called lower and upper solutions. In that sense, they can be considered as founding the lower and upper solutions method. A drawback of Dragoni's approach is that assuming the existence of ordered lower and upper solutions α and β hides the difficulty. In practical problems, there is no clue about finding these functions. This drawback motivated further work to indicate how they can be found. Nagumo [214] proved the existence of at least one solution of the boundary value problem (148) assuming that there exist two functions α and β ∈ C 1 ([a, b]) such that α < β on and assumed that f : E → R is continuous together with the partial derivatives ∂ f ∂u , ∂ f ∂v , and satisfies what today is called the Nagumo condition: for some positive continuous function ϕ : R + → R + that satisfies In fact, Nagumo generalized earlier ideas of Bernstein [215] to control the derivative of the solution. Nagumo's condition, as well as Bernstein's, control the growth of f as a function of the derivative. Although Nagumo's condition is not optimal, it was quite successful due to the fact that it is both general and simple.
Following the upper and lower solution technique due to Chaplygin [216], the Russian school studied the monotone method in some systematic way. Babkin [217] considered the problem (147) and generalized the result of Picard [212] by replacing the decrease of f by a one sided Lipschitz condition and using the approximations (α n ) n obtained as solution of −α n (t) + λ(t)α n (t) = − f (t, α n−1 ) + λ(t)α n−1 (t), α n (a) = α n (b) = 0.
An important step is due to Kantorovich [218] who developed an absolute formulation of the method. This idea was further developed by several others.
Independently of the Russian school, Courant and Hilbert [219] proposed a monotone iterative scheme which was similar to Babkin [217]. The main problem was to find appropriate conditions on the function f to apply the method. Shampine [220] introduced a one-sided Lipschitz condition This condition unifies Courant and Hilbert and other approaches in the literature. Amann [221] generalized the one-sided Lipschitz condition by imposing a Holder condition on f .
The study of the monotone iterative methods for nonlinearities depending on the derivative was initiated by Gendzhoyan [222] who considered the problem u (t) + f (t, u, u ) = 0, u(a) = 0, u(b) = 0.
In the same way he considered the periodic and and Neumann problems.
In all these results discussed above, the usual order (α ≤ β) for the upper and lower solution is considered. In the opposite case, the situation is quite different. The first existence results in the presence of a non-well-ordered lower and the upper solution seems to be due to Amann et al. [225]. A first contribution to the monotone iterative method in case the lower and upper solutions in the reverse order, i.e., α ≥ β was made by Omari and Trombetta [226].

Future Scope
In the following, we list a few problems that researchers may find interesting to address in the near future.
1. ADM is an iterative method that depends on Adomian's polynomial. One may change the polynomial for better accuracy (see Remark 3). 2. One may consider the Equations (82) and (83) as well as (91) and (92) to test a newly developed numerical technique like fuzzy transform finite difference method, nonstandard finite difference methods, iterative methods, etc.
3. One may consider Equation (90) and analyze it to explore problems arising in epitaxial growth problems. This class of problems is recently proposed by Carlos et al. [14] and the references therein. 4. Coupled system of nonlinear SBVPs are not fully explored, so there are various possibilities related to the system of nonlinear SBVPs. 5. Several theoretical as well as numerical works described in Section 2.1.1 may be explored for the case when the nonlinear term is derivative dependent.

Conclusions
The literature related to nonlinear singular DE is huge and difficult to cover. We have made sincere efforts to cover the most important result related to nonlinear singular DE. Though recently there has been less activity in the area of analytical study, several research articles are published related to computation, the reason being that singular BVPs pose a challenge to researchers to compute solutions near the point of singularity. We successfully achieved our goal by writing this survey in a very systematic way. In addition, we present the role of boundary conditions that can help researchers to impose appropriate boundary conditions. When the nonlinear term is of the form f (x, y), we see that lots of theoretical and numerical works are done in comparison to the case when the nonlinear term is of the form f (x, y, py ). We found relatively less theoretical as well as numerical work corresponding to the case when the nonlinear term is of the form f (x, y, py ). We conclude that lots of theoretical and numerical investigations are still pending on a class of nonlinear singular DE with f (x, y, py ) as a nonlinear term, coupled systems when the nonlinear term is of the form f (x, y) or f (x, y, py ). Our review article may help researchers and may save their time, effort, and cost in identifying valid research problems in this area.

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