J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 Journal of the Nigerian Society of Physical Sciences Original Research Perturbed Collocation Method For Solving Singular Multi-order Fractional Differential Equations of Lane-Emden Type O. A. Uwaherena,∗, A. F. Adebisib, O. A. Taiwoa aDepartment of Mathematics, University of Ilorin, Ilorin, Nigeria bDepartment of Mathematics, Osun State University, Oshogbo, Nigeria Abstract In this work, a general class of multi-order fractional differential equations of Lane-Emden type is considered. Here, an assumed approximate solution is substituted into a slightly perturbed form of the general class and the resulting equation is collocated at equally spaced interior points to give a system of linear algebraic equations which are then solved by suitable computer software; Maple 18. Keywords: Perturbed Collocation Method, Legendre polynomial, Singular multi-order differential equations, Lane-Emden equations Article History : Received: 17 March 2020 Received in revised form: 15 May 2020 Accepted for publication: 16 May 2020 Published: 01 August 2020 c©2020 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction Perturbation Collocation Method is a technique of numeri- cal approximation. It has been described as a very useful tool for solving differential and integro-differential equations of dif- ferent kinds [1, 2, 3]. Fractional order differential equation is a special kind of differential equations with fractional order derivative say α, which is a non-integer. Many phenomena in science and engineering resulting in fractional order differen- tial equations have been successfully represented using math- ematical models in the field of mechanics, elasticity, science, education to mention but a few. Most often, the equations aris- ing from the mathematical modelling of fractional order dif- ferential equations are difficult to solve analytically. This is because many of them do not have solutions in closed form ∗Corresponding author tel. no: +23480xxxx572 Email address: uwaheren.ao@unilorin.edu.ng (O. A. Uwaheren ) and hence seeking an approximate solution by numerical treat- ment becomes helpful. To achieve this, a lot of numerical meth- ods have been advanced. [4, 5] proposed Variational Iteration Method (VIM) and used the method to solve multi-order frac- tional integro-differential equations and singular IVPs of Lane- Emden type respectively. Their results showed that the method is an accurate one that yields the exact solution within a few iterations. Their study, however, noted that in some cases, the method may require more calculations which will add some dif- ficulties. [6] solved nonlinear multi-order fractional differential equations using Legendre Pseudo-Spectral Method (LPSM). The results reveal that the method is effective as a small number of shifted Legendre polynomials were needed to obtain a satis- factory result. [7] presented a numerical solution method for solving fractional differential equations using Bernstein poly- nomials. Other methods used by many researchers are Cheby- shev Wavelet method [8], Chebyshev Polynomial Method [9, 10] Collocation Method [11], Pade Approximate Method [12], 141 Uwaheren et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 142 etc. [12] described the Lane-Emden differential equation as a singular initial value problem relating to second-order differen- tial equations which have been used to model several phenom- ena in mathematical physics and astrophysics. According to [13], Lane-Emden type of differential equation is used to describe a variety of phenomena in physics and astro- physics, including isothermal gas spheres and thermionic cur- rents and determining their numerical solutions is very chal- lenging because of the singularity behaviours at the point of origin. The study also solved some fractional order differential equations of Lane-Emden type using the collocation method and admitted that a collocation method is a suitable tool for solving a class of Lane-Emden type of differential equations. [14] worked on the numerical solution of Lane-Emden type equation using multilayer perceptron neural network method. The study aimed at producing an optimal solution of Lane- Emden equation with less computation using multilayer percep- tron artificial neural network technique. The results obtained show that the technique can develop into an effective approach for solving Lane-Emden type problems with less computational time and memory space. In this study, we are concerned with multi-order fractional dif- ferential equations of Lane-Emden type and their solutions by perturbation collocation method. To transform the problems to a system of linear algebraic equations, the collocation method algorithms as given by [13] was implemented. Other researchers that have solved different problems of Lane- Emden type using different methods include [15, 16, 17 ]. The general form of the singular multi-order fractional differ- ential equation is given as n∑ i=0 Piy (i)(x) + Dα (y(x)) + k xα−β Dβ (y(x)) + k xα−2 (y(x)) = f (x) (1) subject to conditions y(0) = A, y′(0) = B (2) where A, B are constants, Pi, (i = 0, 1, 2, ..n) are arbitrary constants to be determined, n is the highest integer order deriva- tive, α and β are fractional order derivatives, and α > β. 2. Definition of Relevant Terms 2.1. Assumed Approximate Solution An Assumed Approximate Solution also called Trial Solu- tion refers to a predetermined polynomial taken to represent the solution of the problem been solved. 2.2. Maple 18 Is a mathematical software that is capable of Performing mathematical computation, for example, solving an equation symbolically or numerically, or construct Maple objects etc de- pending on the command. 2.3. Fractional Differential Equation A differential equation is called a fractional differential equa- tion if it contains at least one fractional order derivatives, Dα of the unknown function y(x). The general form of a fractional differential equation is given as Dαy(x) = f (x, y(x)) (3) subject to the conditions: Dαy(0) = ωk, k = 0, 1, ..., n where Dα is the fractional order derivative in the Caputo sense and α is a non-integer value, n = dαe, called the ceiling α, α > 0 is the highest order of the equation. The operator Dα is defined as Dα∗ f (x) = 1 Γ(m −α) ∫ x 0 (x − t)m−α−1 dm dtm f (t)dt (4) in the Caputo sense for function f (x) with (m−1) absolute con- tinuous derivatives and Riemann-liouville defined it as c D β x f (x) = 1 Γ(n −β) dn d xn ∫ x c (x − t)β−1 f (t)dt (5) The advantage of the Caputo’s fractional derivative is that one can specify the initial conditions of the fractional differential equation in classical form. i.e. D(k)y(k) = bk, k = 0, 1, 2, · · · , m − 1 (6) The generalized factorial form of non-integer order derivatives in Euler’s Gamma function f (x) = xm is given as dα d xα xm = Γ(m + 1) Γ(m − n + 1) xn−m (7) Equation (7) is a property of fractional derivatives given in Euler gamma functions which makes the definition more suit- able and compatible for application. 2.4. Multi-order Fractional differential Equations of Lane-Emden type A fractional multi-order differential equation of Lane-Emden type is of the form equation (1). However, when the highest or- der derivative of the equation is fractional, then equation (1) becomes: Dα (y(t)) + k tα−β Dβ (y(t)) + k tα−2 (y(t)) = f (t) (8) subject to conditions y(0) = A, y′(0) = B (9) where A, B, α and β are as described earlier. 142 Uwaheren et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 143 3. Methodology Consider a general class of multi-order fractional differen- tial equation of Lane-Emden type given in equation (1) together with the conditions in equation (2). In order to solve equations (1) and (2), an assumed approximate solution of the form: yN (t) = N∑ j=0 a j L j(t) (10) is taken. Where L j(t) is Legendre polynomial and N is the de- gree of the assumed approximant. Thus equation(10) is substituted into a slightly perturbed equa- tion (1) to give n∑ i=0 Pi  N∑ j=0 a j L ( j) j (t)  + Dα  N∑ j=0 a j L j(t)  + k tα−β Dβ  N∑ j=0 a j L j(t)  + ktα−2  N∑ j=0 a j L j(t)  = f (t) + Hn(t) (11) where Hn(t) = n∑ (r=1) τr T(N−r−1)(t) (12) is called the perturbation term, n is as stated earlier and T (t) is Chebyshev polynomial. Putting equation (12) into (11), gives n∑ i=0 Pi  N∑ j=0 a j L ( j) j (t)  + Dα  N∑ j=0 a j L j(t)  + k tα−β Dβ  N∑ j=0 a j L j(t)  + ktα−2  N∑ j=0 a j L j(t)  = f (t) + n∑ (r=1) τr T(N−r−1)(t) (13) subject to conditions n∑ k=1 y(k)N (0) = φk; k = 1, 2, ...n (14) Where τr (r = 1(1)n) are the free tau parameters to be de- termined, T(N−r−1)(t) are the Chebyshev polynomials and a j are the unknown constants also to be determined. Applying equations (7) on (13) accordingly gives n∑ i=0 Pi  N∑ j=0 a j L ( j) j (t)  + ktα−β N∑ j=0 a j Γ( j + 1) Γ( j −α + 1) t j−α + k tα−2 N∑ j=0 Γ( j + 1) Γ( j −α + 1) t j−α + N∑ j=0 a j Γ( j + 1) Γ( j −β + 1) t j−β = f (t) + n∑ (r=1) τr T(N−r−1)(t) (15) Equation (15) is further simplified and then collocated at equally spaced interior points, t = te on [a,b]; te = a + (b−a)i (N) ; i = 1,2,...,N to obtain a system of linear algebraic equations, in- cluding those obtained from the use of initial conditions. The system of linear algebraic equations are solved using Gaussian elimination method (using a computer package; Maple 18) to obtain the unknown constants. The constants obtained are then substituted back into the assumed approximate solution to give the required approximate solution. 4. Numerical Examples Example 1 Consider the multi-order fractional differential equation Dα (y(x)) + k x(α−β) Dβ (y(x)) + k x(α−2) (y(x)) = f (x) (16) where = [ 6x ( Γ(4 −β) + kΓ(4 −α) Γ(4 −α)Γ(4 −β) + x2 6 ) −2 ( Γ(3 −β) + kΓ(3 −α) Γ(3 −α)Γ(3 −β) + x2 2 )] x2−α (17) for k = 1, α = 32 and β = 1 2 Subject to the conditions y(0) = y′(0) = 0 0 ≤ x ≤ 1 The exact solution is y(x) = x3 − x2 Equation (16) is perturbed as D 3 2 (y(x)) + 1 x( 3 2 − 1 2 ) D 1 2 (y(x)) + 1 x( 3 2 −2) (y(x)) = f (x) + n∑ r=1 τr T(N−r−1)(x) (18) Taking an approximate solution for N = 4 y4 = 4∑ j=0 y(x) = 4∑ j=0 a j L j(x) = a0 + a1 (2 x − 1) + a2 ( 6 x2 − 6 x + 1 ) + a3 ( 20 x3 − 30 x2 + 12 x − 1 ) + a4 ( 70 x4 − 140 x3 + 90 x2 − 20 x2 + 1 ) (19) Substituting equations (17) and (19) into (18), we have Dα  4∑ j=0 y(x)  + 1x(α−β) Dβ  4∑ j=0 y(x)  + 1x(α−2)  4∑ j=0 y(x)  = 143 Uwaheren et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 144[ 6x( Γ(4 −β) + Γ(4 −α) Γ(4 −α)Γ(4 −β) + x2 6 − 2( Γ(3 −β) + Γ(3 −α) Γ(3 −α)Γ(3 −β) + x2 2 ) ] x2−α + n∑ r=1 τr T(N−r−1)(x) (20) Further simplifications with α = 32 and β = 1 2 gives Dα(a0 +a1(2x−1)+a2(6x 2 −6x+1)+a3(20x 3 −30x2 +12x−1) + a4(70x 4 − 140x3 + 90x2 − 20x + 1)) + 1 x Dα(a0 + a1(2x − 1) + a2(6x 2 −6x + 1) + a3(20x 3 −30x2 + 12x−1) + a4(70x 4 −140x3 + 90x2 − 20x + 1)) + 1 x −1 2 (a0 + a1(2x − 1) + a2(6x 2 − 6x + 1) + a3(20x 3 −30x2 + 12x−1) + a4(70x 4 −140x3 + 90x2−20x + 1)) − ( 6 x ( 28 15 √ π + 1/6 x2 ) − 20 3 √ π − x2 ) √ x −τ1 ( 8 x2 − 8 x + 1 ) −τ2 (2 x − 1) = 0 (21) Applying the operator Dα on equation (21) and after some simplifications gives −1.7724538509 x9/2 + 1.7724538509 x7/2− 56 x5/2 5 + 20 x3/2 3 +124.0717696 x11/2a4+35.44907702 x 9/2a3+10.63472311 x 7/2a2 +3.544907702 x5/2a1+1.7724538509 x 3/2a0−3.544907702 τ2 x 2 −14.17963081 x3τ1−248.1435391 x 9/2a4−53.17361553 x 7/2a3 +1311.520847 x7/2a4−10.63472311 x 5/2a2+245.2694462 x 5/2a3 −1603.449077 x5/2a4−1.7724538509 x 3/2a1+41.77245385 x 3/2a2 −201.7724539 x3/2a3 +601.7724539 x 3/2a4 +6 √ xa1−18 √ xa2 + 36 √ xa3 −60 √ xa4 + 1.7724538509 τ2 x + 14.17963081 x 2τ1 − 1.7724538509 xτ1 = 0 (22) Equation (22) is collocated at N − 2 equal spaced points in the interval [0, 1] and the remaining two equations are obtained from the use of the conditions attached. This gives 7 system of linear equations. Solving the system of equations, we get y4(x) = −1.48141262×10 −10 + 3.2×10−9 x−1.000000039x2 + 1.00000004x3 − 5.469888349 × 10−9 x4 (23) When the same example was solved α = 52 , β = 3 2 following the same methodology,the approximate solution obtained is y4(x) = 5.00 × 10 −11 − 3 × 10−10 x − 1.397902893 x2 + 1.664907851 x3 − 0.1566880882 x4 (24) Example 2 Consider the multi-order fractional differential equation Dα (y(x)) + k x(α−β) Dβ (y(x)) + k x(α−2) (y(x)) = [ 2 ( Γ(3 −β) + kΓ(3 −α) Γ(3 −α)Γ(3 −β) + x2 2 ) − 6x ( Γ(4 −β) + kΓ(4 −α) Γ(4 −α)Γ(4 −β) + x2 6 )] x2−α (25) for k = 1, α = 32 and β = 1 2 subject to conditions y(0) = y′(0) = 0 0 < x ≤ 1 with exact solution y(x) = −x3 + x2 Taking an assumed approximate solution for N = 4, and com- puting the α = 32 and β = 1 2 -derivatives of equation (25) and other necessary simplifications and substitution, gives 35.44907702 x9/2a3+10.63472311 x 7/2a2+3.544907702 x 5/2a1 + 1.7724538509 x3/2a0−3.544907702 τ2 x 2 −14.17963081 x3τ1 −248.1435391 x9/2a4−53.17361553 x 7/2a3+1311.520847 x 7/2a4 −10.63472311 x5/2a2+245.2694462 x 5/2a3−1603.449077 x 5/2a4 −1.7724538509 x3/2a1+41.77245385 x 3/2a2−201.7721531 x 3/2a3 +1.7724538509 x9/2+124.0717696 x11/2a4+601.7724539 x 3/2a4 + 6 √ xa1 −18 √ xa2 + 36 √ xa3 −60 √ xa4 + 1.7724538509 τ2 x + 14.17963081 x2τ1 − 1.7724538509 xτ1 − 1.7724538509 x 7/2 − 20 x3/2 + 51 x2/3 = 0 (26) Equation (26) is collocated at N − 2 equal spaced points in the interval [0, 1] and the remaining two equations are obtained from the use of the conditions attached. This gives 7 system of linear equations. Solving the system of equations, we get y4(x) = 3.77824 × 10 −12 − 1.× 10−10 x + 1.000000109x2 − 0.9999984717x3 − 0.1415135523e − 5x4 (27) When the same example was solved α = 52 , β = 3 2 following the same method,the approximate solution obtained is y4(x) = 6.2 × 10 −9 − 1.449037318 x2 − 8.450067824 x3 + 5.265645764 x4 (28) Example 3 Consider the multi-order fractional differential equation Dα (y(x)) + 2 x Dβ (y(x)) + xy(x) = x2 + x3 + x4 + 96 √ x √ π + 3 √ π √ x ( 6 + 1 √ x ) (29) subject to conditions y(0) = y′(0) = 0 0 < x ≤ 1 with exact solution y(x) = x + x2 + x3 Taking an assumed approximate solution for N = 4, and com- puting the α = 52 and β = 3 2 -derivatives of equation (29) and other necessary simplifications and substitution, gives 70 x5a4 + 20 x 4a3 − 140 x 4a4 + 6 x 3a2 − 30 x 3a3 + 1490 x 3a4 + 2 x2a1 − 6 x 2a2 + 252 x 2a3 − 1700 x 2a4 + xa0 − xa1 + 37 xa2 144 Uwaheren et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 145 −181 xa3 +541 xa4 +4 a1−12 a2 +24 a3−40 a4 = x 4 +8 x3τ1 + x 3 − 8 x2τ1 + 2 x 2τ2 + 12 x 2 + xτ1 − xτ2 + xτ3 + x (30) Equation (30) is collocated at N − 2 equal spaced points in the interval [0, 1] and the remaining two equations are obtained from the use of the conditions attached. This gives 8 system of linear equations. Solving the system of equations, we get y4(x) = 3.2749743×10 (−10) + 1.6×10(−8) x + 1.109848402x2 + 1.000005386x3 − 0.2049775180e − 5x4 (31) When the same example was solved α = 32 , β = 1 2 following the same method,the approximate solution obtained is y4(x) = −4.9638971×10 (−10) + 2.6×10(−8) x + 1.48521236x2 + 1.000023102x3 − 0.346404728e − 5x4 (32) Example 4 Consider the multi-order fractional differential equation Dα (y(x)) + 8 x Dβ (y(x)) + xy(x) = x5 − x4 + 44x2 −30x(33) subject to conditions y(0) = y′(0) = 0 0 < x ≤ 1 with exact solution y(x) = x4 − x3 Taking an assumed approximate solution for N = 4, and com- puting the α = 52 and β = 3 2 -derivatives of equation (33) and other necessary simplifications and substitution, gives − 140 x5a4 + 6 x 4a2 − 30 x 4a3 + 370 x 4a4 + 2 x 3a1 − 6 x 3a2 + 72 x3a3 + 1800 x 3a4 + x 2a0 − x 2a1 + 13 x 2a2 + 419 x 2a3 − 3179 x2a4 + 2 xa1 + 9 xa2 + 70 x 6a4 + 20 x 5a3 − 468 xa3 + 142 xa4 + 16 a1 − 48 a2 + 96 a3 − 160 a4 = x6 −x5 + 8 x3τ1 + 44 x 3 −8 x2τ1 + 2 x 2τ2 −30 x 2 + xτ1 −xτ2 (34) Equation (34) is collocated at N − 2 equal spaced points in the interval [0, 1] and the remaining two equations are obtained from the use of the conditions attached. This gives 8 system of linear equations. Solving the system of equations, we obtained y4(x) = −8.× 10 (−11) + 1.× 10(−9) x − 2.857225303x2 + 4.713873487x3 − .9046244956x4 (35) When the same example was solved α = 32 , β = 1 2 following the same method,the approximate solution obtained is y4(x) = −1.×10 (−10) x−1.000000001x3+1.000000000x4(36) 5. Conclusion In this study, the proposed method was used to solve multi- order fractional differential equations successfully. Four exam- ples were solved and the results of the numerical solutions pre- sented in Tables 1, 2, 3 and 4. When α = 32 and β = 1 2 the pro- posed method produced accurate approximate solution curve Figure 1. Graphical Representation of Error in Table 1 Figure 2. Graphical Representation of Error in Table 2 Figure 3. Graphical Representation of Error in Table 3 which is exactly over the curve of the exact solution for exam- ples 1, and 2. But, the results are not as satisfactory at α = 52 and β = 32 as seen in Figures 1 and 2. However, for examples 3 and 4, there is a change of results. Here, the results are better at α = 52 and β = 3 2 than at α = 3 2 and β = 1 2 as shown in Figure 3 and 4. This is likely because k in these examples is higher than 1 which is 2 and 8 in examples 3 and 4 respectively. The value 145 Uwaheren et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 146 Table 1. Error of Results for Example 1 x α = 52 α = 3 2 β = 32 β = 1 2 Exact Appx Error Appx Error 0.0 0.00000 0.00000 5.0000e-11 -0.000000 5.9232e-13 0.1 0.00900 -0.01232 3.3298e-03 0.008999 2.6668e-10 0.2 0.03200 -0.04284 1.0848e-02 0.031999 9.6994e-10 0.3 0.06300 -0.08213 1.9128e-02 0.062999 1.9604e-09 0.4 0.09600 -0.12112 2.5122e-02 0.095999 3.0781e-09 0.5 0.12500 -0.15115 2.6155e-02 0.124999 4.1532e-09 0.6 0.14400 -0.16393 1.9932e-02 0.144000 5.0056e-09 0.7 0.14700 -0.15153 4.5298e-03 0.146999 5.4456e-09 0.8 0.12800 -0.10640 2.1596e-02 0.127999 5.2733e-09 0.9 0.08100 -0.02139 5.9613e-02 0.080999 4.2787e-09 1.0 0.00000 0.11032 1.1032e-01 0.000000 2.2421e-09 Table 2. Error of Results for Example 2 x α = 52 α = 3 2 β = 32 β = 1 2 Exact Appx Error Appx Error 0.0 0.000000 0.000000 6.2000e-10 0.000000 3.7782e-12 0.1 0.009000 -0.022414 3.1414e-02 0.0090000 2.4706e-09 0.2 0.032000 -0.117137 1.4914e-01 0.0320000 1.4306e-08 0.3 0.063000 -0.315914 3.7891e-01 0.0630000 3.9585e-08 0.4 0.096000 -0.637850 7.3385e-01 0.0960000 7.8988e-08 0.5 0.125000 -1.089415 1.2144e+00 0.125000 1.2980e-07 0.6 0.144000 -1.664440 1.8084e+00 0.1440001 1.8590e-07 0.7 0.147000 -2.344120 2.4911e+00 0.1470002 2.3778e-07 0.8 0.128000 -3.097010 3.2250e+00 0.1280002 2.7253e-07 0.9 0.081000 -3.879029 3.9600e+00 0.0810002 2.7386e-07 1.0 0.000000 -4.633459 4.6335e+00 0.0000002 2.2207e-07 Table 3. Error of Results for Example 3 x α = 32 α = 5 2 β = 12 β = 3 2 Exact Appx Error Appx Error 0.0 0.0000000000 -0.0000000005 4.9639e-10 0.0000000000 3.4707e-11 0.1 0.0110000000 0.0158521485 4.8521e-03 0.0110000001 8.6920e-11 0.2 0.0480000000 0.0674086784 1.9409e-02 0.0480000002 2.1491e-10 0.3 0.1170000000 0.1606697155 4.3670e-02 0.1170000004 3.8136e-10 0.4 0.2240000000 0.3016353775 7.7635e-02 0.2240000006 5.5642e-10 0.5 0.3750000000 0.4963057740 1.2131e-01 0.3750000007 7.1780e-10 0.6 0.5760000000 0.7506810062 1.7468e-01 0.5760000009 8.5069e-10 0.7 0.8330000000 1.0707611670 2.3776e-01 0.8330000009 9.4779e-10 0.8 1.1520000000 1.4625463410 3.1055e-01 1.1520000010 1.0093e-09 0.9 1.5390000000 1.9320366030 3.9304e-01 1.5390000010 1.0431e-09 1.0 2.0000000000 2.4852320250 4.8523e-01 2.0000000010 1.0642e-09 146 Uwaheren et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 147 Table 4. Error of Results for Example 4 x α = 32 α = 5 2 β = 12 β = 3 2 Exact Appx Error Appx Error 0.0 0.0000000000 -0.0000000001 8.0000e-11 0.0000000000 0.0000e+00 0.1 -0.0009000000 -0.0239488420 2.3049e-02 -0.0009000000 1.1000e-11 0.2 -0.0064000000 -0.0780254233 7.1625e-02 -0.0064000000 2.8000e-11 0.3 -0.0189000000 -0.1372031514 1.1830e-01 -0.0189000001 5.7000e-11 0.4 -0.0384000000 -0.1786265321 1.4023e-01 -0.0384000001 1.0400e-10 0.5 -0.0625000000 -0.1816111705 1.1911e-01 -0.0625000002 1.7500e-10 0.6 -0.0864000000 -0.1276437696 4.1244e-02 -0.0864000003 2.7600e-10 0.7 -0.1029000000 -0.0003821324 1.0252e-01 -0.1029000004 4.1300e-10 0.8 -0.1024000000 0.2143448386 3.1674e-01 -0.1024000006 5.9200e-10 0.9 -0.0729000000 0.5285371464 6.0144e-01 -0.0729000008 8.1900e-10 1.0 0.0000000000 0.9520236894 9.5202e-01 -0.0000000010 1.1000e-09 Figure 4. Graphical Representation of Error in Table 4 of the exact solution and the source functions count. From nu- merical results, it can be seen that the proposed method is an accurate estimate for the class of differential equations consid- ered. Acknowledgments The authors wish to sincerely appreciate Mr Akorede M. T, for his immense contributions in the typesetting of this work. We are also grateful to the referees and editor for their invalu- able comments and suggestions. References [1] S. Abbasbandy. “Application of He’s homotopy perturbation method for Laplace transform”, Chaos Solution & Fractals 30 (2006) 1206. [2] M. Dehghan & F. Shakeri. “Solution of an integro-differential equation arising in Oscillating Magnetic Fields using He’s Homotopy Perturbation Method”, PIER 78 (2008) 361. [3] Z. Odibat & S. Momani. “Modified homotopy perturbation method: ap- plication to quadratic riccati differential equation of fractional order”, Chaos, Solution and Fractals 36 (2008) 167. [4] I. K. Omar. “Variational iteration method for solving multi fractional in- tegrodifferential equations”, Iraqi Journal of Science 55 (2014) 1086. [5] A. Yildirim & T. Qzis. “Solution of singular IVPs of Lane-Emden type by variational iteration method”, Nonlinear Analysis: Theory, Method & Applications 70 (2009) 2480. doi.org/10.1016/j.na.2008.03.012 [5] Y. Yang. “Solving a nonlinear multi-order fractional differential equation using Legendre pseudo-spectral method”, Journal of Applied Mathemat- ics 4 (2013) 113. [6] L. Jianping, L. Xia, X. Liyong & S. Yufa. “Solution of fractional differen- tial equations by Bernstein polynomials method”, Journal of Convergence Information Technolgy 7 (2012) 820. [7] A. Barzkar, M. Oshagh, P. Assari & M. Mehrpouya. “Numeri- cal solution of the nonlinear Fredholm integral equation and Fred- holm integro-differential equation of second kind using Cheby- shev wavelets”, World Applied Sciences Journal 18 (2012) 1774. doi.org/10.5829/idosi.wasj.2012.18.12.920. [8] Z. Avazzadeh & M. Heydari. “Chebyshev polynomials for solving two dimensional linear and nonlinear integral equations of the second kind”. Computational and Applied Mathematics 31 (2012) 127. [9] H. Cerdik-Yaslan & Z. Akyuz-Dascioglu. “Chebyshev polynomial solu- tion of nonlinear Fredholm-Volterra integro-differential equation”, Jour- nal of Arts Science 5 (2006) 89. [10] O. A. Uwaheren & O. A. Taiwo. “A collocation technique based on an orthogonal polynomial for solving multi-fractional order integro differen- tial equations”, Journal Mathematical Society of Nigeria (ABACUS) 43 (2016) 215. [11] M. Yigider. “The numerical method for solving differential equations of Lane-Emden type by Pade approximation”, Discrete Dynamics in Nature and Society (2011) 1. doi.org/10.1155/2011/479396 [12] S. M. Mohammed & S. Norazak. “Numerical study of fractional differ- ential equations of Lane-Emden type by method of collocation”, Applied Mathematics 3 (2012) 851. doi.org/10.4236/am.2012.38126. [13] A. Verma & M. Kumar. “Numerical solution of Lane-Emden type equation using multilayer perceptron neural network method”, Interna- tional Journal of Applied and Computational Mathematics 5 (2019) 141. link.springer.com [14] A. H. Bhrawy & S. A. Alofi. “A Jacob Gauss collocation method for solv- ing nonlinear Lane-Emden type equation”, Communications in Nonlinear Science and Numerical Simulation 17 (201) 62. [15] K. Parand, M. Dehgham, A. Rezaeia & S. Ghaderi. “An applica- 147 Uwaheren et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 141–148 148 tion algorithm for the solution of the nonlinear Lane-Emden type equations arising in astrophysics using Hermite functions colloca- tion method”, Computer Physics Communications 181 (2010) 1096. doi.org/10.1016/j.cpc.2010.02.018 [16] W. Al-Hayani, L. Alzubaidy & A. Entesar. “Solutions of singular IVPs of Lane-Emden type by homotopy analysis method with genetic algo- rithm”, Applied Mathematics and Information Sciences 11 (2019) 1. doi.org/10.18576/amis/paper 148