A Successive Linearization Method Approach to Solve Lane-Emden Type of Equations

Abstract

We propose a new application of the successive linearization method for solving singular initial and boundary value problems of Lane-Emden type. To demonstrate the reliability of the proposed method, a comparison is made with results from existing methods in the literature and with exact analytical solutions. It was found that the method is easy to implement, yields accurate results, and performs better than some numerical methods.

1. Introduction

In this work, we investigate the application of a novel approach of solving a class of nonlinear singular initial value and boundary value problems in the second-order ordinary differential equations. In recent years, these problems have attracted the attention of many researchers because of their useful applications in astronomy, mathematical biology, physics, and other areas of science and engineering. One of the equations describing this type of equations is the Lane-Emden type equation formulated as 𝑦^′′(𝑥)+2𝑥𝑦^′(𝑥)+𝑦^𝑛=0, (1.1) with initial conditions 𝑦(0)=1,𝑦^′(0)=0, (1.2) where 𝑛∈[0,5] is a constant parameter. This equation is very useful in astrophysics in the study of polytropic models and stellar structures [1, 2]. For the special case when 𝑛=0,1,5, exact analytical solutions were obtained by Chandrasekhar [1]. For all other values of 𝑛, approximate analytical methods and numerical methods are used to approximate the solution of the Lane-Emden equation. Analytical approaches that have recently been applied in solving the Lane-Emden equations include the Adomian decomposition method [3, 4], differential transformation method [5], homotopy perturbation method [6], homotopy analysis method [7, 8], power series expansions [9–13], and variational iteration method [14, 15]. Generally, when all the above cited analytical approaches are used to solve Lane-Emden equation, a truncated power series solution of the true solution is obtained. This solution converges rapidly in a very small region (0<𝑥<1). For 𝑥>1, convergence is very slow and the solutions are inaccurate even when using a large number of terms [16]. Convergence acceleration methods such as Padé approximations may be used to improve the convergence of the resulting series or to enlarge their domain of convergence. The homotopy analysis method [8, 17] has a unique advantage over the other analytic approximation methods because it has a convergence controlling parameter that can be adjusted to improve the region of convergence of the resulting series.

An important physical parameter associated with the Lane-Emden function is the location of its first positive real zero. The first zero of 𝑦(𝑥) is defined as the smallest positive value 𝑥₀ for which 𝑦(𝑥₀)=0. This value is important because it gives the radius of a polytropic star. The analytic approaches on their own are not very useful in solving for 𝑥₀ because their region of convergence is usually less than 𝑥₀. Recently, there has been a surge in the number of numerical methods that have been proposed to find solutions of the Lane-Emden equations. Recent numerical methods that have been proposed include the Legendre Tau method [18], Sinc-collocation method [19], the Lagrangian approach [20], the successive linearization method [21], and a numerical method based on radial basis functions [22]. Accurate results for the Lane-Emden equation have previously been reported in [23] where the Runge-Kutta routine with self-adapting step was used to generate seven digit tables of Lane-Emden functions. These tables are now widely used as a benchmark for testing the accuracy of new methods for solving the Lane-Emden equations. More recently, numerical values giving the first zero of the Lane-Emden equation for various values of 𝑛 with at least ten decimal places of accuracy were tabulated in [24] who use a numerical perturbation approach and in [25] who use the Chebyshev spectral method.

Another important class of singular ODE's of Lane-Emden type is the Frank-Kamenetskii [26] boundary value problem 𝑦^′′(𝑥)+𝑘𝑥𝑦^′(𝑥)+𝛿𝑒^𝑦=0, (1.3) which is solved subject to the boundary conditions 𝑦^′(0)=0,𝑦(1)=0. (1.4) This equation is used to model the temperature distribution in a vessel before a thermal explosion takes place. A substantial amount of research has been done on this type of problem using both analytical and numerical treatments. Closed form solutions of this problem were reported in [26, 27] for the case when 𝑘=0 and 1. Analytical treatment of (1.3) leading to series solutions was carried out in [28–30] using the Adomian decomposition method and in [31] using a power series approach based on symmetry reduction. Asaithambi [32] used a shooting method based approach that uses automatic differentiation to obtain a Taylor series expansion of the solution. Numerical approaches based on B-spline functions (see e.g, [28, 33, 34]) have been widely used to develop solutions of singular boundary value problems of the type (1.3). Kumar and Gupta [35] present a survey and review of papers of spline solution of singular boundary value problems of the type (1.3).

In the present paper, we present a very simple and robust method that gives very accurate solutions to the nonlinear initial and boundary value problems governed by (1.1) and (1.3), respectively. The method extends the approach of [21] who attempted to solve the Lane-Emden equation (1.1) for the case when 𝑛=2,3,4 using the successive linearization method (SLM). The successive linearization method is a new numerical iterative approach that has been recently reported and successfully utilized in solving boundary value problems [36–39] occurring in boundary layer flows in the application of fluid mechanics.

The SLM approach is based on transforming an ordinary nonlinear differential equation into an iterative scheme made up of linear equations which are then solved using numerical approaches. In [21], the SLM was used to obtain solutions of (1.1) that were accurate only up to seven decimal places. In this study, a modified version of the SLM is presented which converges much faster and gives solutions that are accurate to more than fifteen decimal places. We also extend the application of the present method to singular boundary value problems governed by the Frank-Kamenetskii equation. The accuracy of the present method it tested against results obtained using other methods was found in, literature.

2. Numerical Solution

In this section, we present the method used to solve the governing equations (1.1)–(1.4).

2.1. Solution of Lane-Emden Equation

To solve the Lane-Emden equation (1.1), it is convenient to recast the problem from being an initial value problem to a boundary value problem by considering only the domain 𝑥∈[0,𝛼], where 𝛼 is the first zero of 𝑦(𝑥). In most practical applications of the Lane-Emden equation (1.1), the goal is to integrate the governing equation from 0 to 𝛼. Since 𝛼 is an unknown parameter, we rescale the problem by setting 𝑥=𝛼𝑡. (2.1) Substituting (2.1) into (1.1) and simplifying gives 𝑦_𝑡𝑡+2𝑡𝑦_𝑡+𝛼²𝑦^𝑛=0,𝑦(0)=1,𝑦(1)=0,𝑦_𝑡(0)=0, (2.2) which is a nonlinear eigenvalue problem with 𝛼 as the eigenvalue. To solve (2.2), we use the successive linearization method (SLM) (see e.g., [36–39]). The SLM approach is based on transforming an ordinary nonlinear differential equation into an iterative scheme made up of linear differential equations which are then solved using analytical or numerical methods wherever possible. In applying the SLM on (2.2), we set 𝑦(𝑡)=𝑦_𝑖(𝑡)+^𝑖−1_𝑚=0𝑦_𝑚(𝑡),𝑖=1,2,3,…, (2.3) where 𝑦_𝑖 are unknown functions that are obtained by iteratively solving the linearized version of the governing equations assuming that 𝑦_𝑖(0≤𝑚≤𝑖−1) are known from previous iterations. The algorithm starts with an initial approximation 𝑦₀(𝑡) which is chosen to satisfy the boundary conditions in (2.2). A suitable initial guess that satisfies all the boundary conditions in (2.2) of this example is 𝑦₀(𝑡)=1−𝑡². (2.4) We observe that (2.2) is a second-order differential equation but has three boundary conditions. The third boundary condition (𝑦_𝑡(0)=0) is used as an extra condition that is required to solve for the eigenvalue 𝛼. Since 𝛼 is an unknown parameter, we also apply the SLM approach in solving for 𝛼 by setting 𝛼=𝛼_𝑖+^𝑖−1_𝑚=0𝛼_𝑚. (2.5) An initial approximation for 𝛼 that was found to give accurate results for all 𝑛∈(0,5) was 𝛼²₀4=63^𝑛. (2.6) This value of 𝛼₀ was obtained by substituting the initial approximation 𝑦₀ given by (2.4) in (2.2) and solving the resulting equation at 𝑥=1/2, the middle of the domain for 𝑡.

Substituting the expansions (2.3) and (2.5) into (2.2) and neglecting nonlinear terms in 𝑦_𝑖 and 𝛼_𝑖 give 𝑦̈_𝑖+2𝑡𝑦̇_𝑖+𝑎_1,𝑖−1𝑦_𝑖+𝑎_2,𝑖−1𝛼_𝑖=𝑟_𝑖−1, (2.7) subject to 𝑦_𝑖(0)=0,𝑦_𝑖(1)=0,𝑦̇_𝑖(0)=0, (2.8) where 𝑎_1,𝑖−1=𝑛^𝑖−1_𝑚=0𝛼_𝑚^2𝑛−1_𝑚=0𝑦_𝑚^𝑛−1,𝑎_2,𝑖−1=2^𝑖−1_𝑚=0𝛼_𝑚^𝑛−1_𝑚=0𝑦_𝑚^𝑛,𝑟_𝑖−1=−^𝑖−1_𝑚=0𝑦̈_𝑚+2𝑡^𝑖−1_𝑚=0𝑦̇_𝑚+^𝑖−1_𝑚=0𝛼_𝑚^2𝑖−1_𝑚=0𝑦_𝑚^𝑛, (2.9) where the dots denote differentiation with respect to the scaled variable 𝑡. To solve (2.7), we use the Chebyshev collocation method in which case the functions 𝑦_𝑖 are approximated by Chebyshev interpolating polynomials in terms of their values at the collocation points. Before applying the spectral method, it is convenient to transform the physical region [0,1] of problem 11 into the domain [−1,1]. This can be achieved by using the mapping 𝑡=𝑧+12,−1≤𝑧≤1. (2.10) We use the Gauss-Lobatto collocation points [40, 41] to define the Chebyshev nodes in [−1,1] as 𝑧_𝑗=cos𝜋𝑗𝑁,𝑗=0,1,…𝑁, (2.11) where (𝑁+1) is the number of collocation points. The derivatives of 𝑦_𝑖 at the collocation points can be written as 𝑦̇_𝑖𝑧_𝑗=^𝑁_𝑘=0𝐷_𝑗𝑘𝑦_𝑖𝑧_𝑘=𝐷𝑌_𝑖,𝑦̈_𝑖𝑧_𝑗=^𝑁_𝑘=0𝐷²_𝑗𝑘𝑦_𝑖𝑧_𝑘=𝐷²𝑌_𝑖, (2.12) where 𝑌_𝑖=[𝑦_𝑖(𝑧₀),𝑦_𝑖(𝑧₁),…,𝑦_𝑖(𝑧_𝑁−1),𝑦_𝑖(𝑧_𝑁)]^𝑇, 𝑇 is transpose, and 𝐷 is the Chebyshev derivative matrix of size (𝑁+1)×(𝑁+1) whose entries are defined [40, 41] as 𝐷_𝑗𝑘=𝑐_𝑗𝑐_𝑘(−1)^𝑗+𝑘𝑧_𝑗−𝑧_𝑘𝐷𝑗≠𝑘;𝑗,𝑘=0,1,…,𝑁,_𝑘𝑘𝑧=−^𝑘21−𝑧²_𝑘𝐷𝑘=1,2,…,𝑁−1,₀₀=2𝑁²+16=−𝐷_𝑁𝑁, (2.13) with 𝑐_𝑘=2,𝑘=0,𝑁,1,−1≤𝑘≤𝑁−1. (2.14) Applying the Chebyshev spectral collocation method in (2.7)-(2.8) gives 𝐴_𝑖−1𝑌_𝑖+𝐵_𝑖−1𝛼_𝑖=𝑅_𝑖−1, (2.15) subject to 𝑦_𝑖𝑧_𝑛=0,𝑦_𝑖𝑧₀=0,_𝑁𝑘=0𝐷_𝑁𝑘𝑦_𝑖𝑧_𝑘=0. (2.16)(2.17) The matrices in (2.15) are defined as 𝑅_𝑖−1=𝑟_𝑖−1𝑧₀,𝑟_𝑖−1𝑧₁,…,𝑟_𝑖−1𝑧_𝑁−1,𝑟_𝑖−1𝑧_𝑁^𝑇,𝐴_𝑖−1=4𝐷²+4𝐭_𝑑𝐷+𝐚_1,𝑖−1,𝐚_1,𝑖−1=𝑎_1,𝑖−1𝑡₀0⋯00𝑎_1,𝑖−1𝑡₁⋯0⋮⋮⋮⋮00⋯𝑎_1,𝑖−1𝑡_𝑁,4𝐭_𝑑=4𝑡₀040⋯0𝑡₁4⋯0⋮⋮⋮⋮00⋯𝑡_𝑁,𝐵_𝑖−1=𝐚_2,𝑖−1,with𝐚_2,𝑖−1=𝑎_2,𝑖−1𝑡₀𝑎_2,𝑖−1𝑡₁⋮𝑎_2,𝑖−1𝑡_𝑁. (2.18) We note that the system (2.15) consists of (𝑁+1) unknowns in 𝑌_𝑖 and an additional unknown 𝛼_𝑖. To solve for all (𝑁+2) unknowns, we increase the dimension of the system (2.15) by introducing an additional row on which we impose the Neumann boundary condition (2.17). The resulting system takes the form 𝐴_𝑖−1𝐵_𝑖−1𝐷_𝑁,0⋯𝐷_𝑁,𝑁0𝑌_𝑖𝛼_𝑖=𝑅_𝑖−10. (2.19)

The boundary conditions (2.16) can be imposed on system (2.19) by deleting the first and (𝑁+1)st rows of the coefficient matrix and vectors and deleting the first and (𝑁+1)st columns of the coefficient matrix. Thus, starting from the initial approximations 𝑌₀ and 𝛼₀, the solutions 𝑌_𝑖 and 𝛼_𝑖(𝑖=1,2,3,…) can be determined from 𝑌_𝑖𝛼_𝑖=𝐴_𝑖−1𝐵_𝑖−1𝐷_𝑁,1⋯𝐷_{𝑁,𝑁−1}0⁻¹𝑅_𝑖−10. (2.20)

2.2. Solution of the Frank-Kamenetskii Equation

The Frank-Kamenetskii equation admits exact analytical solutions when 𝑘=1 (see e.g., [27, 42]). When 𝛿<2, it has two solutions. The problem has no solution for 𝛿>2 and when 𝛿=2, it has a unique solution. The solutions are given by [42] 𝑦=log16𝑒^𝑐12𝛿+𝑒^𝑐1𝑥²²,𝑦=log16𝑒^𝑐21+2𝛿𝑒^𝑐2𝑥²², (2.21) where the constants 𝑐₁ and 𝑐₂ are given by 𝑐₁=log2(4−𝛿)±42(2−𝛿),𝑐₂=log4−𝛿±22(2−𝛿)2𝛿². (2.22)

To solve (1.3), it is convenient to introduce the boundary condition at 𝑥=0 as 𝑦(0)=𝛽, (2.23) where 𝛽 is a constant. The governing equation (1.3) is then solved subject to the Dirichlet boundary conditions 𝑦(0)=𝛽 and 𝑦(1)=0 using the SLM approach as described in the previous subsection using the Neumann boundary condition 𝑦^′(0)=0 as an extra condition for determining the unknown 𝛽. To this end, we look for a solution of the form 𝑦(𝑥)=𝑦_𝑖(𝑥)+^𝑖−1_𝑚=0𝑦_𝑚(𝑥),𝛽=𝛽_𝑖+^𝑖−1_𝑚=0𝛽_𝑚, (2.24) where 𝑦_𝑖 and 𝛽_𝑖 are obtained iteratively by solving the linearized equations that result from substituting (2.24) into the governing equation (1.3). Substituting (2.24) into (1.3)-(1.4) and neglecting nonlinear terms in 𝑦_𝑖 and 𝛽_𝑖 give 𝑦_𝑖^′′+𝑘𝑥𝑦^′_𝑖+𝑏_1,𝑖−1𝑦_𝑖=𝑠_𝑖−1, (2.25) subject to 𝑦_𝑖(0)−𝛽_𝑖=𝑝_𝑖−1,𝑦_𝑖(1)=0,𝑦^′_𝑖(0)=0, (2.26) where 𝑏_1,𝑖−1=𝛿exp^𝑖−1_𝑚=0𝑦_𝑚,𝑝_𝑖−1=^𝑖−1_𝑚=0𝛽_𝑚−^𝑖−1_𝑚=0𝑦_𝑚(0),𝑠_𝑖−1=−^𝑖−1_𝑚=0𝑦_𝑚^′′+𝑘𝑥^𝑖−1_𝑚=0𝑦^′_𝑚+𝛿exp^𝑖−1_𝑚=0𝑦_𝑚. (2.27) A suitable initial guess that satisfies the boundary conditions (1.4) and (2.23) is 𝑦₀(𝑥)=𝛽₀1−𝑥², (2.28) where 𝛽₀ is the initial approximation for 𝛽 which is given by 𝛿exp3𝛽₀4=2𝛽₀(1+𝑘). (2.29) Equation (2.29) is obtained by substituting (2.28) into the governing equation (1.3) and setting 𝑥=1/2. The solutions for 𝛽₀ can be obtained using nonlinear equation solvers in scientific computing software such as the fsolve routine that is available in MAPLE or MATLAB.

Equations (2.25)-(2.26) are solved using the Chebyshev spectral collocation method. The details of the implementation of this approach are similar to those given for the Lane-Emden equation in the previous section. The discretized equation system to be solved is 𝐶_𝑖−1𝑌_𝑖=𝑆_𝑖−1,_𝑁𝑘=0𝐷_𝑁𝑘𝑦_𝑖𝑧_𝑘=0,𝑦_𝑖𝑧₀𝑦=0,_𝑖𝑧_𝑛−𝛽_𝑖=𝑝_𝑖−1, (2.30) where 𝑆_𝑖−1=𝑠_𝑖−1𝑧₀,𝑠_𝑖−1𝑧₁,…,𝑠_𝑖−1𝑧_𝑁^𝑇,𝐶_𝑖−1=4𝐷²+2𝑘𝐱_𝑑𝐷+𝐛_1,𝑖−1,𝐛_1,𝑖−1=𝑏_1,𝑖−1𝑥₀0⋯00𝑏_1,𝑖−1𝑥₁⋯0⋮⋮⋮⋮00⋯𝑏_1,𝑖−1𝑥_𝑁,2𝑘𝐱_𝑑1=2𝑘𝑥₀010⋯0𝑥₁1⋯0⋮⋮⋮⋮00⋯𝑥_𝑁. (2.31) The equation system can be written as the matrix equation 0𝐶10⋯000_𝑖−1⋮00⋯01−1𝐷_𝑁,0𝐷_𝑁,1⋯𝐷_{𝑁,𝑁−1}𝐷_𝑁,𝑁0𝑦_𝑖𝑧₀𝑦_𝑖𝑧₁⋮𝑦_𝑖𝑧_𝑁−1𝑦_𝑖𝑧_𝑁𝛽_𝑖=0𝑠_𝑖−1𝑧₁⋮𝑠_𝑖−1𝑧_𝑁−1𝑝_𝑖−10, (2.32) where the boundary conditions have been imposed on the first and last rows of 𝐶_𝑖−1 and 𝑆_𝑖−1 and the derivative boundary condition has been added in the (𝑁+2) row. Thus, starting from 𝑦₀ and 𝛽₀ given by (2.28) and (2.29), the subsequent solutions for 𝑦_𝑖,𝛽_𝑖(𝑖=1,2,3,…) can be obtained by solving the matrix system (2.32).

3. Results

In this section, we present the approximate solutions of the governing equations (1.1)–(1.4) using the successive linearization method (SLM). In order to assess the performance and reliability of the present method of solution, the results are tabulated in Tables 1–6 and compared with accurate results from the literature and exact analytical solutions for the cases where analytical solutions are available.

Table 1 

Comparison of the first zero (𝛼) of the Lane-Emden equation at different iterations for different values of 𝑛 against numerical results given by Boyd [25] and Seidov [24].

Table 2 

Comparison of the first zero (𝛼) of the Lane-Emden equation at different iterations for different values of 𝑛 against numerical results given by Boyd [25].

Table 3 

Comparison of the value of 𝑑𝑦/𝑑𝑡(1), for the Lane-Emden equation at different iterations for different values of 𝑛 against numerical results given by Boyd [25].

Table 4 

SLM results for 𝑦(0) corresponding to the two roots of (1.3) when 𝑘=1, 𝛿=1, and 𝑁=50.

Table 5 

Maximum absolute errors for solution of (1.3) (1st root).

Table 6 

Maximum absolute errors for solution of (1.3) (2nd root).

In the case of the Lane-Emden equation (1.1), exact analytical solutions are only available for the case when 𝑛=0,1,5. For any other values of 𝑛, numerical methods are used to integrate the equation. The main difficulty in analyzing the Lane-Emden equation is the singularity behaviour occurring at 𝑥=0. Most researchers revisit the Lane-Emden equation repeatedly to test the viability of their new numerical methods against it. Accurate numerical results to seven decimal places are tabulated in Horedt [23, 43]. More accurate results including more decimal places have recently been reported in [24, 25]. In this work, we compare the results of the SLM approximate results to at least twelve decimal places for the dimensionless radius and mass of the polytropic stars with the results of [24, 25]. In Table 1, we present the SLM solution of the Lane-Emden equation (1.1) for the first zero 𝛼. The first zero is the smallest root 𝛼 such that 𝑦(𝛼)=0 and gives the dimensionless radius of the polytropic star [1]. The present results, at different iterations, are compared with the numerical results of [24, 25] to fifteen decimal places. The results presented in Table 1 we generated using 𝑁=40 collocation points for 𝑛=1,2 and 𝑁=70 for 𝑛=3,4. It can be seen from the table that the present SLM approach gives very accurate results which rapidly converge to the numerical result. Accuracy to fifteen decimal places is obtained after only four iterations for 𝑛=1, five iterations for 𝑛=2, six iterations for 𝑛=3, and seven iterations for 𝑛=4. We remark that the original SLM solution of the Lane-Emden equation reported in [21] gave results which were only accurate to seven decimal places even after increasing the number of iterations. This shows that the present modification of the original SLM approach results in a more robust method of solution. Similar results for nonintegral values of 𝑛 are given in Table 2.

In Table 3, we give a comparison of the values of 𝑦^′(1) which is used in the definition of the dimensionless mass of the polytropic star given by −𝛼²𝑦^′(𝛼) against the numerical results of [25]. Again, it can be seen that the present SLM results rapidly converge to the numerical results of [25].

Figure 1 shows the Lane-Emden equation solution curves for different indices obtained using eight iterations of the SLM solution series. These solution curves are qualitatively similar to those obtained using other numerical methods (see e.g., [22]).

Figure 1 

Graph of the SLM solutions of the Lane-Emden equation (1.1) for 𝑛=0.5,1,1.5,2,2.5,3,3.5, and 4 using 𝑖=8 iterations.

The Frank-Kamenetskii equation (1.3) was solved using the SLM iteration scheme given by the matrix equation (2.32), and the present results were compared against the exact analytical results given by (2.21)-(2.22) for the case when 𝑘=1 and 𝛿=1. The initial approximations required to start the SLM algorithm were generated using (2.28)-(2.29). For 𝑘=1 and 𝛿, (2.29) generates two solutions for 𝛽₀ which then result in two SLM solutions. Table 4 gives a comparison between the SLM approximate results to twenty decimal places against the exact analytical solutions. It can be seen from the table that for the first root, convergence to twenty decimal places is achieved after only three iterations. For the second root, convergence to twelve decimal places of accuracy is achieved after only five iterations and convergence to twenty decimal places is achieved after six iterations. From this observation, it can be seen that the present algorithm is very robust.

In Table 5, we give results for the maximum absolute errors for (1.3) when 𝑘=1 and 𝛿=1 for different collocation points (grid points) at different iterations. The present results are compared with the recent B-Spline numerical results of [33]. It is worth noting that the B-Spline approach of [33] was meant to be an improvement on other numerical approaches that have previously been used to solve (1.3). For example, the results reported in [33] are more accurate than those reported in [34] where a three point finite difference approach was used to solve the same problem. Table 5 indicates that the present results are far more superior than the B-Spline results of [33] when the same number of grid points are used. We remark that results for only the first root were reported in [33] whereas our approach gives two solutions of (1.3). It can also be seen from Table 5 that the convergence of the method improves with an increase in the number of iterations.

Table 6 gives the results for the second root. Again, it can be seen from this table that the present method gives very accurate results which rapidly converge to the exact solution. Increasing 𝑁 and the number of iterations improves the accuracy of the method.

In Figure 2, we plot the solution profile 𝑦(𝑥) for different values of the Frank-Kamenetskii parameter 𝛿 after 6 iterations. It can be seen from the figure that there is good agreement between our approximate results and the exact analytical solutions.

Figure 2 

Comparison between the analytical solution (in diamonds) and the SLM approximate solutions for different values of 𝛿 when 𝑘=1.

4. Conclusion

In this paper, we presented a new application of the successive linearization method (SLM) in solving Lane-Emden equations that model polytropic models of arbitrary index 𝑛∈(0,5). The method is also tested on a singular boundary value Frank-Kamenetskii problem. The governing nonlinear equations were transformed to eigenvalue problems and subsequently solved using the SLM approach. A comparison was made between exact analytical solutions, numerical results from the literature and the present approximate solutions. The comparison indicates that the present method is robust and gives very accurate results and performs better that other numerical methods that have previously been applied to solve the same problems. The results presented in this paper can easily be extended to other initial and boundary value problems which are difficult to solve using other numerical methods.