Various boundary value problems associated to EF equation have been treated along time by transforming them into initial value problems. Furthermore, they have been solved with different classical numerical methods of various accuracy.

The IVP in physics as well as in other branches of science, have been solved by various shooting methods. These methods have become obsolete over time for many reasons.

We appreciate that solving boundary problems with (their) specific numerical methods is important in order to obtain better outgoings.

Consider the Emden-Fowler equation

supplied with the same boundary condition as above, i.e.,

When $s>1,$ the solution of this equation decays as $x^{-q}$, $q=3/\left(2\left(s-1\right)\right)$, as $x$ approaches $\mathrm{\infty }$ (see [4] and the papers quoted there).

The existence of the solution to such a problem is substantiated in the paper [1]. In [4] the author finds very accurate a priori bounds for the slope $u_0^{\prime }:=u^{\prime }\left(0\right).$ He exploits the integral properties of EF equation and produces the double inequality

Thus, the collocation spectral methods, in conventional form or more recent Chebfun form, more elegant, easier to implement and more accurate are well suited for solving the problems at hand.

In the case of the problems examined in this article, the Laguerre polynomials, multiplied by the weight function $exp(-bx)$ correctly represent the continuous functions on the real semi-axis and satisfy the boundary condition at infinity (the second condition in (2)).

Conventional spectral methods have been exploited almost to exhaustion solving physics problems by J P Boyd [6].

The main aim of this study is to show that the LGRC method solves problems of type (1)–(2) fairly accurate such that the slope in origin satisfies bounds like (3). Moreover we infer on the convergence of this algorithm.

Delkhosh and Parand have recently published several papers on the Emden-Fowler problem from which we only quote [9]. They confirm the efficiency and accuracy of the spectral methods using the generalized fractional order of the Chebyshev orthogonal functions. The purpose of our work is only to eliminate the numerical noise from the value of the derivative of solution in origin. The extension of our study to non-integer values of $s$ is in progress.

For the particular case $s:=3/2$ in [4] the author deduces the following important and useful integral identity

for any natural $n.$

Some clarifications are necessary for the convergence of the improper integral [5].

The solutions of problems (1)–(2) are continuous indefinitely differentiable functions with compact support on the real semi-axis, the integrand of integral I has the same property.

Under these conditions, the general theory of improper integrals ensures convergence.

For $n:=2$ and $n:=4$ we will show that the solutions obtained by LGRC satisfy the above integral identities with a quite good approximation. The final bounds for $u_0^{\prime }$ are in this case

In the paper [6] the author gives very interesting historical aspects related to the numerical solution of the TF problem, mentioning several Nobel Prize laureates who considered it. He lists a set of values for $u_0^{\prime }$, not all of which belong to the range (3). Its solution is based on rational Chebyshev series.

With respect to the LGRC method we refer to the monograph [5] which offers all the ingredients to implement the method.

Actually we are looking for a solution of the form

where

and $x_1=0$ and $x_2,\mathrm{\dots },x_N$ are the roots of $L_{N-1}\left(x\right)$, the Laguerre polynomial of degree $N-1,$ indexed in increasing order of magnitude. The interval $[0,\mathrm{\infty })$ can be mapped to itself by the change of variables $x:=b\tilde{x},$ where $b$ is a positive real number called scaling factor. Thus LGRC method contains a free parameter which can be exploited to optimize the accuracy of differencing process and to tune the numerical outcomes.

In order to introduce the boundary conditions (2) we operate the change of variables

If we denote by ${𝒟}_{N-2}^{\left(2\right)}$ the second order Laguerre differentiation matrix in physical space, LGRC cast the problem (1)–(2) with $s:=3/2$ into the nonlinear algebraic system

where the vector $X_{N-2}$ contains the Gauss-Radau nodes $x_2,\mathrm{\dots },x_{N-1}$ and the vector $U$ the corresponding unknowns (see Weidemann and Reddy [3]). We have to observe that in formulating the system (7) we use the MATLAB notations (the under-dot in (7) signifies the pointwise product of the two areas in MATLAB).

We have solved the system (7) by MATLAB routine fsolve and have found the following value

which strictly belong to the interval (5). In order to find the slope at origin we simply multiply the vector solution with the first order derivative matrix.

In the left panel of Figure 1 we display the solution to our problem and in the right panel we show how the coefficients of LGRC solution decrease. Roughly speaking we cannot hope to an accuracy much better than ${10}^{-4}.$ (see Boyd [2], Convergence of Spectral Methods). Actually, the way in which the coefficients decrease is a correct indication of convergence of solution.

We have to mention at this point that in order to get the coefficients of solution we have relied on the Laguerre polynomial transform between physical and coefficient spaces from [5].

With respect to the right panel in Figure 1, we have to observe that this technique was introduced by Boyd in his huge monograph [2]. It is not a precise proof but a very valuable trick. We have successfully and frequently used this argument.

The Newton’s method in solving the system (7) is convergent as it is visible from Figure 2.

With the value of $u_0^{\prime }$ displayed above we have validated numerically the integral relations (4). Thus, for $n:=2$ the integral equality is correct to order ${10}^{-3}$ and when $n:=4$ the accuracy decreases to order ${10}^{-2}.$

We now consider the EF problem with $s:=4$ in the equation (1) and a perfect similar treatment as above has been applied. For this value of exponent $s$ the bounds in (3) are the following

Our numerical experiments with scaling factor $b$ in interval $[25,100]$ and resolution $N$ around $300$ provided numerical values of the derivative in origin strictly in the interval (8).

The only cons of the method is the arbitrary choosing of the resolution $N$.

In order to solve nonlinear BVP we can choose between Chebyshev or Fourier collocation methods, when the domain is finite, or Laguerre, sinc or Hermite polinomials based collocation, whenever the domain is unbounded. The choosing depends also on the shape of nonliniarity as well as on the packages of implementation at hand. We mention that the literature about these methods is now well established so we consider that more details are not necessary.

In order to validate problems (1)–(2) we have alternatively tried to use Chebyshev collocation in its conventional form (see [7]) or in form of Chebfun (see [8]) with two ingredients, namely domain truncation or the mapping of half line into the canonical interval $[-1,1].$ None of these methods gave acceptable results.

It is clear that the success of the used LGRC method is due to the existence of the weighting factor $exp(-x/2)$ in the definition of the LGR interpolant (6). The method handles correctly both singularities, at origin and at infinity, is robust with respect to the resolution and scaling, and is easy implementable. Up to our knowledge we consider that this method is, so far, the most accurate in solving such nonlinear and double singular problems as boundary value ones.

The author is indebted to both referees for opportune observations that led to the improvement of the work.