There exists a vast literature on domain truncation type methods, the most used methods for problems defined on unbounded domains, but considerable less rigorous estimations of the domain truncation errors. Domain truncation was coupled with shooting, finite differences and Galerkin schemes. Galerkin or collocation methods based on the Laguerre functions, and the Chebyshev polynomials with various mappings, are the other two alternatives (see for instance the paper of Boyd et al. [11]).

Variants of shooting method along with domain truncation have been usually used to solve numerically two-point boundary value problems on unbounded domains, i.e., on the real line or on the half line. For high order problems, the method is based on the use of certain asymptotic initial conditions together with the compound matrix method. Some iterative and non iterative transformations, which reduce the problems to solutions of a sequence of Cauchy’s related ones, have also been considered (see for instance the paper of Fazio [14]). Along the last decades, using all these techniques, a lot of numerical as well as important analytic results have been established.

In some of our recent papers [15] and [16] we successfully used Laguerre collocation (LC for short) in order to avoid the well known drawback of domain truncation. In this paper we use additionally another variant of LC, namely the Laguerre Gauss

Radau collocation (LGRC for short) introduced by Shen et al. in [31] or Benacchio and Bonaventura [8]. This latter variant has the advantage of a discrete transform.

Thus, the main aim of this paper is to solve accurately higher order boundary value problems, including eigenvalues ones, by a collocation method based on Laguerre functions. The method is fairly efficient in imposing the boundary conditions in origin and at infinity, and easily implementable. The boundary conditions at infinity are directly satisfied by Laguerre functions and the conditions at zero are introduced by a removing technique of independent boundary conditions. With this technique the analytical singularities at the origin have been completely resolved.

The free parameter involved in the LC algorithms is adjusted in order to improve the convergence of the Newton type processes which solve the nonlinear algebraic systems generated by LC and LGRC discretizations. A favorable choice of parameter could speed up the process and enlarge the attraction basin of the solution.

Moreover, we additionally make use of solutions provided by LC as a way to validate and extend some analytical results of existence, uniqueness and asymptotic behavior concerning the solutions to second and third orders nonlinear boundary value problems on the half line. With respect to some challenging BVPs from physics and engineering we show that generally all the analytical results are accurately validated by this pseudospectral method. More important is the fact that we can get a better understanding of the singularities of problems. For instance we detect boundary layer behavior of solutions when the analytical result fails to notice this behavior. However, there are also some exceptions. For instance, in a boundary layer problem, (see Sect. 3.1) our existence results can not be overlapped on some analytic older ones.

All the designed MATLAB codes are fairly similar and simple and do not require large memory space. They work efficiently in CPU times of order of a few seconds or less.

The outline of this paper is as follows: in Section 2 we first briefly review the LC and LGRC algorithms. We have improved our previous numerical results concerning the gas diffusion into an semi-infinite porous medium in the subsection 2.1 (the so called Kidder problem). We confirm that the analytic singularity is a week one. Then we solve two benchmark problems. A closed solution problem which arise in a variety of fields viz. plasmas, colloids, semiconductors, etc. is analyzed in subsection 2.2. The dependence of the accuracy of solutions on the scaling factor is here noticed. An eigenvalue problem attached to Schrödinger equation, with various potentials, is considered in subsection 2.3 . Using the non normality of the left hand side matrix in the pencil (Henrici number and pseudospectrum) we heuristically justify the efficiency of the LC method. In the first subsection of Section 3 we solve a boundary layer problem encountered in physical oceanography. A singular plasma physics problem is analyzed in the subsection 3.2 and in subsection 3.3 we solve a classical ground state problem on the half line. Some remarks on the performances of LC algorithm are provided in the Section 4.

Just in a very few words we recall the main features of the LC method exposed for instance in our booklet [17], Ch. 4. We consider second and third order problems of the form

where $L^{(n)}$ is a linear $n$th order ordinary differential operator, $f$ is an arbitrary (nonlinear) function of $u(x)$ and its $n-1$ derivatives and $\lambda$ is a real (physical) parameter. Actually, in our applications we will consider the cases $n=2$ and 3 . The maps $g_{k}$ are some simple linear combination of the values of $u$ and its first and possibly second derivatives in zero and at infinity. Eigenvalue problems are supposed to belong to this class of problems.

Existence and uniqueness results, as well as results concerning the asymptotic behavior of solutions, for classes of problems belonging to (1) have been reported in the literature, see for instance the papers of Agarwal [4] and O’Regan [29]. A particular class of problems belonging to (1) are the positive radial solutions of some Dirichlet problems (see for instance [24] and [25] to quote but a few).

In order to solve the two-point boundary value problem (1) by LC we represent $u(x)$ by the interpolate $p_{N-1}(x)$ defined by

where $\phi_{j}(x)$ are the Lagrangian interpolating polynomials in barycentric form

The interpolating nodes $x_{j},j=2,\ldots,N$ are the first $N-1$ roots of Laguerre polynomial $L_{N-1}(x)$ of degree $N-1$ indexed in increased order of magnitude. We add the node $x_{1}=0$ in order to facilitate the incorporation of boundary conditions. We recall that $x_{N}=O(N)$ as $N\rightarrow\infty$.

The interval $[0,\infty)$ can be mapped to itself by change of variables $x=b\widetilde{x}$, where $b$ is any positive real number (the scaling factor). The LC method therefore contains a free-parameter. It means that the Laguerre interpolating process is exact for functions of the form

where $p(x)$ is any polynomial of degree $N-1$ or less.
The usual Laguerre interpolation was systematically analyzed in some well known textbooks (see for instance [9]).

However, the most important result remains that provided in [9] which reads

where $I_{N}$ is the usual interpolation operator and the weights $\omega$ are defined as $\omega_{0}:=e^{-x}$ and $\omega_{\tau}:=e^{-(1-\tau)x},0<\tau<1$.

In order to write down the LC equations for (1) we use the Laguerre differentiation matrices from [34] denoted by $\mathbf{D}^{(n)}$, where $n=1,2,3$ is the differentiation order.

Unfortunately, these matrices are fully populated, non symmetric and quite bad conditioned.

With respect to the boundary conditions, the boundary conditions at infinity are directly satisfied by Laguerre functions. In fact we have

for any natural number $m$.
The conditions at zero are introduced by a removing technique of independent boundary conditions introduced in [21]. In a few words, with this technique we consider the boundary conditions on the origin as $l$ linear constraints on the state of our system of $N$ degrees of freedom and remove them since they are slaved to the other $N-l$ degrees which we keep. Thus, we decompose the state of the system in degrees of freedom that we keep and degrees of freedom that we remove. We express the removed state variables as a function of the kept ones, thus preserving their action in the system. Moreover, we are not interested in the evolution of the removed degrees of freedom since we can, at any time, recover them from the kept ones.

The LC casts the problem (1) into the following nonlinear algebraic system

where $\widetilde{\mathbf{U}}$ contains the nodal unknown values $u_{j},j=l+1,\ldots,N$, i.e. the kept degrees of freedom, and the matrices $\widetilde{\mathbf{D}}^{(k)},k=1,\ldots,n$ are square matrices of the dimension $N-l$. The tilde over these matrices means that the first $l$ rows and columns of the entire differentiation matrices are deleted and the boundary conditions have already been enforced. The vector $X$ contains the kept nodes $x_{l+1},\ldots,x_{N}$. The boundary conditions

are automatically fulfilled due to the linearity of maps $g_{k}$.
The Laguerre Gauss Radau quadrature nodes and weights are defined for instance in the monograph of J. Shen et al. [31] p. 243. The set of nodes contains the origin
and in ascending order the zeros of the derivative of the generalized Laguerre polynomial of degree $N$. On these nodes the authors define the differentiation matrix, say again $\mathbf{D}$, and then along with the above weights they introduce the discrete Laguerre transform and its inverse. The (forward) discrete transform will be fairly useful in order to transfer the information (values of solution in nodes) from the physical space into the frequency (phase) space. Thus we get the LGRC expansion coefficients of the solution and can quantify them (see Sect. 2.1). Actually, the LGRC casts the problem (1) into a nonlinear algebraic system fairly similar to (7). Observe that in this case $\mathbf{D}^{(m)}=\mathbf{D}^{m}$.

Solutions to various nonlinear systems (7) have been carried out using the MATLAB built in function fsolve as well as using a discrete transcript of the very efficient and specialized quasilinearization algorithm (QL for short) provided by Mandelzweig and Tabakin in [2]. The QL iterative algorithm is a generalization of the classical Newton-Raphson method. For problems of the class (1), under fairly general conditions on the nonlinearities $F$ and $G$, it provides in a small number of iterations fast convergent (quadratic), monotonic (often) and stable numerical results. These conditions refer to the existence of functional derivatives $f_{u^{(k)}},k=0,\ldots n-1$ such that the mean value theorem is applicable.

For practical purposes, it is important to underline that the zeroth iteration has to be chosen from physical and mathematical considerations. Usually, it is advantageous that this initial iteration would satisfy at least one of the boundary conditions. All these requirements are fulfilled for the examples below.

The scaling factor has been adjusted manually based on the convergence of the algorithms just mentioned.

Throughout this paper the scaling factor involved in the Laguerre collocation methods has been manually adjusted. This has been done with respect to the convergence outcomes provided by the MATLAB solver fsolve (exitflag and the residual) as well as with regard to the residual provided by the quasilinearization algorithm. An a priori determination or an automatic adaptation of this free parameter, for a specific problem, is still an open issue.

In spite of this minor inconvenience the LC algorithms proved to be fairly simple and accurate in solving a set of genuinely nonlinear and, at the same time, singular boundary value problems defined on the half line. Neither an arbitrary truncation of this domain, nor large systems of ordinary differential equations, supposed by compound matrix methods have been necessary. All the nonlinear algebraic systems encountered within this paper have reasonable dimensions (of order $N$ ) and the CPU elapsed time, required to accurately solve them, does not exceed some few seconds on an ordinary machine.

Additionally, the method has offered real clues in order to clarify the dependence of solutions on some physical parameters and with the purpose of accurately resolve some analytic singularities. Comparison with some other numerical and analytical results from literature, shows that both Laguerre collocation methods are highly competitive and can be used as companion tools for various theoretical studies.

Acknowledgement. The author thanks to the anonymous referee for carefully reading the manuscript. She/he pointed out an annoying flaw. The author is also sincerely indebted to the managing editor for her efficient and generous work.

[1] R. P. Agarwal, D. O’Regan, Boundary Value Problems on the Half Line in the Theory of Colloids, Math. Probl. Eng., 8(2002), 143-150.
[2] R. P. Agarwal, D. O’Regan, Infinite Interval Problems Modeling the Flow of a Gas Through a Semi-Infinite Porous Medium, Stud. Appl. Math., 108(2002), 245-257.
[3] R. P. Agarwal, D. O’Regan, Infinite interval problems modelling phenomena which arise in the theory of plasma and electrical potential theory, Stud. Appl. Math., 111(2003), 339-358.
[4] R. P. Agarwal, Singular Boundary Value Problems With Real World Applications, (2008) http://www.dmmm.uniroma1.it/ agostino.prastaro/AGARWAL-LECTURE.pdf Accessed 22 Feb 2015
[5] H. Alici, H. Taşeli,The Laguerre pseudospectral method for the radial SchrĂśdinger equation, Appl. Numer. Math., doi:10.1016/j.apnum.2014.09.001
[6] O. Bayrak, I. Boztosun, Analytical solutions of the HulthĂŠn and the Morse potentials by using the Asymptotic Iteration Method, J. Mol. Structure: THEOCHEM, 802(2007), 17-21.
[7] J. V. Baxley, Existence and Uniqueness for Nonlinear Boundary Value Problems on Infinite Intervals, J. Math. Anal. Appl., 147(1990), 122-133.
[8] T. Benacchio, L. Bonaventura, Absorbing boundary conditions: a spectral collocation approach, Int. J. Numer. Meth. Fluids (2013), DOI: 10.1002/fld. 3768
[9] C. Bernardy, Y. Maday, Spectral Methods, In Ciarlet, P, Lions, JL., (eds.) Handbook of Numerical Analysis, V. 5 (Part 2), North-Holland, 1997.
[10] L. E. Bobisud, Existence of Positive Solutions to some Nonlinear Singular Boundary Value Problems on Finite and Infinite Intervals, J. Math. Anal. Appl., 173(1193), 69-83.
[11] J.P. Boyd, C. Rangan, P.H. Bucksbaum, Pseudospectral methods on a semi-infinite interval with application to the hydrogen atom a comparison of the mapped Fourier sine method with Laguerre series and rational Chebyshev expansions, J. Comput. Phys. 188(2003), 56-74.
[12] C. V. Coffman, Uniqueness of the ground state solution for $\Delta u-u+u^{3}=0$ and a variational characterization of other solutions, Arch. Rational Mech. Anal., 46(1972), 12-95.
[13] M. Countryman, R. Kannan, Nonlinear Boundary Value Problems on Semi-Infinte Intervals, Comput. Math. Appl., 28(1994), 59-75.
[14] R. Fazio, A novel approach to the numerical solution of boundary value problems on infinite intervals, SIAM J. Numer. Anal., 33(1996), 1473-1483.
[15] C. I. Gheorghiu, Laguerre collocation solutions to boundary layer type problems, Numer. Algorithms 64(2013) 385-401.
[16] C. I. Gheorghiu, Pseudospectral solutions to some singular nonlinear BVPs. Applications in nonlinear mechanics, Numer. Algor., 68(2015), 1-14.
[17] C. I. Gheorghiu, Spectral Methods for Non-Standard Eigenvalue Problems. Fluid and Structural Mechanics and Beyond, Springer, Cham Heidelberg New York Dordrecht London, 2014.
[18] C. I. Gheorghiu, J. Rommes, Application of Jacobi-Davidson mathod to accurate analysis of singular hydrodynamic stability problems, Int. J. Numer. Meth. Fluids, 71(2013), 358-369.
[19] C. I. Gheorghiu, M. E. Hochstenbach, B. Plestenjak, J. Rommes, Spectral collocation solutions to multiparameter Mathieu’s system, Appl. Math. Comput., 218(2012), 11990-12000.
[20] R. Hammerling, O. Koch, C. Simon, E. B. Weinmüller, Numerical solution of singular eigenvalue problems for ODEs with a focus on problems posed on semi-infinite intervals, ASC Report No. 8/2010 Vienna University of Technology, Vienna, Austria
[21] J. Hoepffner, Implementation of boundary conditions, (2010) http://www.lmm.jussieu.fr/ hoepffner/research/realizing.pdf Accessed 2 Feb 2015
[22] R. Iacono, J. P. Boyd, The Kidder Equation: $u_{xx}+2xu_{x}/\sqrt{1-\alpha u}=0$., Stud. Appl. Math., 135(2014), 63-85.
[23] V. B. Mandelzweig, F. Tabakin, Quasilinearization approach to nonlinear problems in physics with applications to nonlinear ODEs, Comput. Phys. Commun., 141(2001), 268-281.
[24] K. McLeod, J. Serrin, Uniqueness of positive radial solutions of $\Delta u+f(u)=0$ in $\mathbb{R}^{n}$. Arch. Rational Mech. Anal., 99(1987), 115-145.
[25] K. McLeod, W. C. Troy, F. B. Weissler, Radial Solutions of $\Delta u+f(u)=0$ with Prescribed Numbers of Zeros, J. Diff. Eq., 83(1987), 368-378.
[26] G. R. Ierley, O. G. Ruehr, Analytic and numerical solutions of a nonlinear boundary layer problem, Stud. Appl. Math., 75, (1986), 1-36
[27] G. R. Ierley, Boundary Layers in General Ocean Circulation, Annu. Rev. Fluid Mech., 22(1990), 111-142.
[28] B. Plestenjak, C. I. Gheorghiu, M. E. Hochstenbach, Spectral collocation for multiparameter eigenvalue problems arising from separable boundary value problems, J. Comput. Phys., 298(2015), 585-601.
[29] R. O’Regan, Solvability of Some Singular Boundary Value Problems on the Semi-Infinite Interval, Can. J. Math., 48(1996), 143-158.
[30] A. Roy, The generalized pseudospectral approach to the bound states of the HulthĂŠn and the Yukawa potentials, PRAMANA-J. Phys., 65(2005), 1-15.
[31] J. Shen, T. Tang, L-L. Wang, Spectral Methods. Algorithms, Analysis and Applications, SpringerVerlag Berlin Heidelberg, 2011.
[32] W. C. Troy, Solutions of Nonlinear Boundary Layer Problems Arising in Physical Oceanography, Rocky Mt. J. Math., 21(1991), 813-820.
[33] W. C. Troy, The existence and uniqueness of bound-state solutions of a semi-linear equation, Proc. R. Soc. A., 461(2005), 2941-2963.
[34] J. A. C. Weideman, S. C. Reddy, A MATLAB Differentiation Matrix Suite, ACM Trans. on Math. Software, 26(2000), 465-519.