\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx, amssymb} \usepackage{array} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 331, pp. 1--18.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/331\hfil Solving Thomas-Fermi equation] {A new approach for solving nonlinear Thomas-Fermi equation based on fractional order of rational Bessel functions} \author[K. Parand, A. Ghaderi, H. Yousefi, M. Delkhosh \hfil EJDE-2016/331\hfilneg] {Kourosh Parand, Amin Ghaderi, Hossein Yousefi, Mehdi Delkhosh} \address{Kourosh Parand (corresponding author) \newline Department of Computer Sciences, Shahid Beheshti University, G.C., Tehran, Iran.\newline Department of Computer Sciences, Department of Cognitive Modelling, Institute for Cognitive and Brain Sciences, Shahid Beheshti University, G.C., Tehran, Iran} \email{k\_parand@sbu.ac.ir} \address{Amin Ghaderi \newline Department of Computer Sciences, Shahid Beheshti University, G.C., Tehran, Iran} \email{amin.g.ghaderi@gmail.com} \address{Hossein Yousefi \newline Department of Computer Sciences, Shahid Beheshti University, G.C., Tehran, Iran} \email{hyousefi412@gmail.com} \address{Mehdi Delkhosh \newline Department of Computer Sciences, Shahid Beheshti University, G.C., Tehran, Iran} \email{mehdidelkhosh@yahoo.com} \thanks{Submitted June 16, 2016. Published December 27, 2016.} \subjclass[2010]{34B16, 34B40, 74S25} \keywords{Fractional order of rational Bessel functions; Collocation method; \hfill\break\indent Thomas-Fermi equation; Quasilinearization method; Semi-infinite domain; Nonlinear ODE} \begin{abstract} In this article, we introduce a fractional order of rational Bessel functions collocation method (FRBC) for solving the Thomas-Fermi equation. The problem is defined in the semi-infinite domain and has a singularity at $x = 0$ and its boundary condition occurs at infinity. We solve the problem on the semi-infinite domain without any domain truncation or transformation of the domain of the problem to a finite domain. This approach at first, obtains a sequence of linear differential equations by using the quasilinearization method (QLM), then at each iteration the equation is solves by FRBC method. To illustrate the reliability of this work, we compare the numerical results of the present method with some well-known results, to show that the new method is accurate, efficient and applicable. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction}\label{sec1} Many problems in mathematics, fluid dynamics, quantum mechanics, astrophysics, physics, and engineering are arisen on the infinite or semi-infinite domains. In this section, we have expressed some of the approaches for solving problems which are defined in unbounded domains and a brief history of Thomas-Fermi equation that is defined on the semi-infinite domain. \subsection{Solving problems over unbounded domains} Recently, various approaches have been successfully proposed for solving problems which are arisen on unbounded domains. Such as numerical, analytical and semi-analytical methods. Different numerical methods have been introduced to the problems which is defined in the semi-infinite domain, such as the Finite difference method (FDM) \cite{refdd1,refdd2}, Finite element method (FEM) \cite{refdd2,refdd3}, Meshfree methods \cite{refdd4,refdd5,refdd5_1}, and Spectral methods \cite{refdd6,refdd8}. The study of analytical and semi-analytical solutions of differential equations (DEs) plays an important role in mathematical physics, engineering, and the other sciences. In the past several decades, various methods for obtaining solutions of DEs have been presented, such as the Adomian decomposition method \cite{refdd19,refdd20}, Homotopy perturbation method \cite{refdd21,refdd211}, Variational iteration method \cite{refdd22}, Exp-function method \cite{refdd23} and so on. The spectral approximations for DEs on finite domains have achieved great success and popularity in recent years, but spectral approximations for DEs on infinite domains have only received limited attention. Several spectral methods for treating infinite/semi-infinite domain problems have been utilized by different researchers: (1) There is an effective approach to solve these problems by applying the basis of Sinc, Hermite and rational Christov functions that are orthogonal on the interval $(-\infty,\infty)$ \cite{refdd9,refdd10} and the basis of Laguerre polynomials that are orthogonal on the interval $[0,\infty)$ \cite{refdd11,refdd12}. (2) Another approach for solving such problems which is based on rational approximations. This method transfers polynomials on interval $[\alpha,\beta]$ to functions on interval $[0,\infty)$ by using the algebraic mapping $x\to \frac {\beta x+\alpha L}{x + L}$ that $L>0$ is a scaling/stretching factor \cite{refdd14}. The Jacobi polynomials are a class of classical orthogonal polynomials, also the Gegenbauer polynomials, the Legendre and Chebyshev polynomials, are special cases of these polynomials which have been used in several literatures for solving some problems. In some papers have been provided the collocation method for natural convection heat transfer equations embedded in porous medium, nonlinear differential equation and nonlinear integro-differential equation based on rational Gegenbauer, Legendre, Chebyshev functions \cite{refdd8,refdd11,refdd24,refdd25}. Doha et al have presented Jacobi rational-Gauss collocation method based on Jacobi rational functions and Gauss quadrature integration to solve nonlinear Lane-Emden equation \cite{refdd16}. Isik et al have used Bernstein polynomials to solve high order initial and boundary value problems. Their approximate solution has a better convergence rate than the one found by using the collocation method \cite{refdd15}. (3) Guo \cite{Guo1,Guo2} has applied a method that proceeds by mapping the original problem in an unbounded domain to a problem in a bounded domain, and then using suitable Jacobi polynomials such as the Gegenbauer polynomials to approximate the resulting problems. (4) A further approach consists of replacing the infinite domain with $[-K,K]$ and the semi-infinite domain with $[0,K]$ by choosing $K$ sufficiently large. This method is named domain truncation \cite{reff001,refdd17}. In this investigation, we attempt to introduce a Spectral method based on the fractional order of rational Bessel functions (FRB) to solve Thomas-Fermi on the semi-infinite domain. \subsection{Thomas-Fermi equation} One of the most important nonlinear ordinary differential equations that occurs in semi-infinite interval is Thomas-Fermi equation as follows \cite{refd00,refd01,refd02}: \begin{equation}\label{eqqq01} \frac{d^2y(x)}{dx^2}-\frac{1}{\sqrt{x}}y^{3/2}(x)=0, \quad x\in[0,\infty), \end{equation} where the boundary conditions for this equation are as follows: \begin{equation}\label{eqqq02} y(0)=1,\quad \lim_{x\to\infty}y(x)=0. \end{equation} The Thomas-Fermi equation appears in the problem of determining the effective nuclear charge in heavy atoms, and because of its importance to theoretical physics, computing its solutions has attracted the attention of the Nobel laureates John Slater (chemistry) \cite{refd03} and Richard Feynman (physics) \cite{refd04} and of course Enrico Fermi \cite{refd05}. One measure of the rapidity of the convergence of the procedure is provided by the calculation of the value of the initial slope $y'(0)$ of the Thomas-Fermi potential \cite{refd06}. The problem is useful for calculating form-factors and for obtaining effective potentials which can be used as initial trial potentials in self-consistent field calculations. The initial slope $y'(0)$ is difficult to compute by any means and plays an important role in determining many physical properties of the Thomas-Fermi atom. It determines the energy of a neutral atom in the Thomas-Fermi approximation: \begin{equation}\label{eqqq03} E=\frac{6}{7}\Big(\frac{4\pi}{3}\Big)^{2/3}Z^{7/3}y'(0), \end{equation} where $Z$ is the nuclear charge. For these reasons, the problem has been studied by many researchers and by the different techniques have been solved, that a number of them are as follows: Baker in 1930 \cite{refd07} studied the singularity of this equation and calculated an analytical solution as follows: \[ y(x)=1-Bx+\frac{4}{3}x^{3/2}-\frac{2}{5}Bx^{5/2}+\frac{1}{3}x^3 +\frac{3}{70}B^2x^{7/2}-\frac{2}{15}Bx^4+\dots, \] where $-B$ is the value of the first derivative at the origin that has calculated $y'(0)=-B=-1.588558$. Esposito in 2002 \cite{refd14} reported an original method, due to Majorana, that leads to a semi-analytical series solution of the Thomas-Fermi equation with appropriate boundary conditions in terms of only one quadrature, and proved that the series expansion is uniformly convergent in the interval $[0,1]$, and has calculated $y'(0)=-1.588$. Liao in 2003 \cite{refd16} employed the Homotopy analysis method and gave an explicit analytic solution of the Thomas-Fermi equation and the related recurrence formula of constant coefficients. The corresponding $m$th-order approximation is \[ y(x)=\sum_{k=0}^{m}\sum_{n=1}^{4k+1}\alpha_{k,n}(1+x)^{-n}, \] where $\alpha_{k,n}$ is defined in \cite[Eq. 26]{refd16}. He calculated $y'(0)=-1.58712$. Kobayashi et al in 1955 \cite{refd31} examined the asymptotic solution of obtained by Coulson and March \cite{refd28}, and improved their solution: \[ y(x)=\frac{144}{x^3}\big(1-z+0.6256974977z^2-0.3133861150z^3 +0.1373912767z^4-\dots \big), \] where $z=\frac{F}{x^c}$, and $F=13.27097391$ and $c=0.7720018726$, and calculated $y'(0)=-1.588070972$. Adomian in 1998 \cite{refh1} introduced a standard decomposition method for solving Thomas-Fermi equation. As briefly as follows: \begin{equation} y(x)=c_{1}+c_{2}x+L^{-1}x^{-1/2}\sum_{n=0}^{\infty}A_{n} \end{equation} where $L^{-1}$ denotes a two-fold integration, $A_n$ denotes the Adomian polynomials generated for $y^{3/2}$, and $c_1$, $c_2 $ are constants of integration. The Adomian decomposition method employs the recursive relation. Marinca and Herianu in 2011 \cite{refh4} used a new method to find an analytical approximate solution to Thomas-Fermi equation and called it the Optimal Parametric Iteration Method (OPIM) that this new iteration approach provides us with a convenient way to optimally control the convergence of the approximate solution. This new iteration approach containing a new iteration scheme involves the presence of a finite number of initially unknown parameters, which are optimally determined. In this way, the approximate initial slope is $y'(0)=-1.5880659888022421$. Zhu et al in 2012 \cite{refh5} approximated the original Thomas-Fermi equation by a nonlinear free boundary value problem (FBVP) and applied an iterative method to solve the FBVP. They transformed the FBVP to a nonlinear singular BVP defined on $[0, 1]$ by a change of variables and also employed an adaptive finite element method based on moving mesh to obtain the best approximate solution at each iteration. Best approximation obtained by this method is $y'(0)=-1.58794357$. A simple and more precise solution to the Thomas-Fermi equation is obtained by making use of the famous Ritz Variational method. Oulne in 2011 \cite{refh7} used a new simple Variational solution of the Thomas-Fermi equation which reproduces the numerical solution accurately in a wide range with a correct asymptotic behavior at long distances from the origin and which allows us to calculate with exactness the initial slope. The proposed solution will be developed in power series which have the same form as series solutions that have been obtained previously by Baker \cite{refd07}. In this method, the approximate initial slope is $y'(0)=-1.588071034$. Abbasbandy and Bervillier in 2011 \cite{reff10} compared three methods based respectively on Taylor (Maclaurin) series, Pad\'{e} approximates and conformal mappings. $y'(0) = -1.5880710226113753127189 \pm 7 * 10^{-22}$ was obtained by using the Pad\'{e}-Hankel method. Boyd in 2013 \cite{reff01} applied collocation method based on the rational Chebyshev functions on semi-infinite intervals $TL_n(y; L)$ which $L$ is a user-choosable numerical. Boyd employed Newton-Kantorovich iteration to reduce the nonlinear differential equation to a sequence of linear differential equations, and he has calculated$$y'(0) = -1.5880710226113753127186845$$with $L=64$ and $600$ collocation points. MacLeod in 1992 \cite{reff03} used two differing approximations on Chebyshev polynomial according to behavior Thomas-Fermi function, one for small $x<40$, one for large $x$. In this method, the approximate initial slope is $y'(0)=-1.5880710226$. Parand et al \cite{reff055,reff05,reff07,reff09} proposed collocation method on rational Chebyshev, Hermite polynomials and Sinc functions to solve Thomas-Fermi on semi-infinite interval without truncating it to a finite domain. These methods reduce the solution of this problem to the solution of a system of algebraic equations. Jovanovic et al in 2014 \cite{reff06} solved the Thomas-Fermi equation by applying a spectral method using an exponential basis set in a semi-infinite domain. The goal of the spectral method approach is to find the values of coefficients $a_i$ that best satisfy the equation \begin{equation} y(x)=\sum_{i=1}^{N} a_iR_i ,\quad R_i=e^{-\beta_i x}, \end{equation} where values of $R_i$ are selected in an intuitive way to cover all the possible decay rates. They have reported detailed about the convergence rate of the initial slope $y'(0)$ for an exponential basis set. Liu and Zhu in 2015 \cite{reff08} proposed an iterative method based on the Laguerre pseudospectral approximation which the solution of the Thomas-Fermi equation as the sum of two parts due to its singularity at the origin. One ``singular'' part is a power series expansion. The other ``smooth'' part satisfies a nonlinear two-point boundary value problem. In this method, the approximate initial slope is $y'(0)=-1.588072$. Yao in 2008 \cite{refm12} solved the Thomas-Fermi equation with a kind of analytic technique, named Homotopy analysis method and his answer is $y'(0)=-1.588004950$. Amore et al in 2014 \cite{refm16} obtained highly accurate solutions to the Thomas-Fermi equations for atoms and atoms in very strong magnetic fields. And they apply the Pad\'{e}-Hankel method, numerical integration, power series with Pad\'{e} and Hermite-Pad\'{e} approximates and Chebyshev polynomials. They solved Thomas-Fermi for different $x$ and obtain answers for $y(x)$ and $y'(x)$. Their best answer is $y'(0)=-1.588071022611375312718684509$. Fernandez in 2011 \cite{refm17} showed that a simple and straightforward rational approximation to the Thomas-Fermi equation provides the slope at the origin with unprecedented accuracy and that Pad\'{e} approximates of relatively low order are far more accurate than more elaborate approaches proposed recently by other authors. He calculated $y(x)$ for different values of $x$ and compare their method with Chebyshev and numerical method, and calculated $y'(0)=-1.588071022611375313$. Epele et al in 1999 \cite{refm18} used Pad\'{e} approximate approach to solving Thomas-Fermi equation. They have calculated $y'(0)=-1.5881$. Khan and Xu in 2007 \cite{refm19} used an analytic technique, namely the Homotopy analysis method (HAM). Their best answer for $y'(0)$ was $-1.586494973$ when they selected [30,30] for Homotopy-Pad\'{e} approximations. The rest of this paper is arranged as follows: Section \ref{sec2} introduces a novel the fractional order of rational Bessel functions (FRB). Section \ref{sec3} describes a brief formulation of quasilinearization method (QLM) introduced in \cite{refff05}. In section \ref{sec4} at first, by utilizing QLM over Thomas-Fermi equation a sequence of linear differential equations is obtained and then at each iteration the fractional order of rational Bessel functions collocation method (FRBC) is used for solving the linear differential equations. We in section \ref{sec5} compared our solutions with some well-known results, comparisons show that the present solutions are highly accurate, we also describe our results via tables and figures. Finally, we give a brief conclusion in section \ref{sec6}. \section{Fractional order of rational Bessel functions $(FRB)$}\label{sec2} The Bessel functions arise in many problems in physics possessing cylindrical symmetry, such as the vibrations of circular drumheads and the radial modes in optical fibers. Bessel functions are usually defined as a particular solution of a linear differential equation of the second order which known as Bessel's equation. Bessel functions first defined by the Daniel Bernoulli on heavy chains (1738) and then generalized by Friedrich Bessel. More general Bessel functions were studied by Leonhard Euler in (1781) and in his study of the vibrating membrane in (1764) \cite{reffff01,reffff02}. \subsection{Definition of Bessel polynomials} The Bessel differential equation of order $n\in\mathbb{R}$ is \begin{equation}\label{eqqq04} x^2\frac{d^2y(x)}{dx^2}+x\frac{dy(x)}{dx}+(x^2-n^2) y(x)=0, \quad x\in(-\infty,\infty). \end{equation} One of the solutions of equation \eqref{eqqq04} by applying the method of Frobenius as follows \cite{reffff03}: \begin{equation}\label{eqqq05} J_{n}(x)=\sum_{r=0}^{ \infty}\frac{(-1)^r}{r!(n+r)!}(\frac{x}{2})^{2r+n}, \end{equation} where series \eqref{eqqq05} is convergent for all $x\in(-\infty,\infty)$. Bessel functions and polynomials are used to solve a number of problems in physics, engineering, mathematics, and etc., such as Blasius equation, Lane-Emden equations, integro-differential equations of the fractional order, unsteady gas equation, systems of linear Volterra integral equations, high-order linear complex differential equations in circular domains, systems of high-order linear Fredholm integro-differential equations, etc. \cite{reffff04,reffff07,reffff08,reffff09,reffff10,reffff11,reffff13,reffff1302}. Bessel polynomials have been introduced as follows \cite{reffff12}: \begin{equation}\label{eqqq06} B_{n}(x)=\sum_{r=0}^{[\frac{N-n}{2}]}\frac{(-1)^r}{r!(n+r)!} (\frac{x}{2})^{2r+n},\quad x\in[0,1]. \end{equation} where $n\in\mathbb{N}$, and $N$ is the number of the basis of Bessel polynomials. Let $\Gamma=\{x:0\leq x \leq 1 \}$ and $L^{2}_{w}(\Gamma)=\{v :\Gamma \to \mathbb{R}| v$ is measurable and $\| v \|_{w} < \infty \}$, where $$ \| v \|_{w}=\Big(\int^{1}_{0}|v(x)|^{2}w(x)dx\Big)^{1/2}, $$ with $w(x)=1$, is the norm induced by the inner product of the space $L^{2}_{w}(\Gamma)$ as follows: $$ \langle v(x),u(x)\rangle_{w}=\int^{1}_{0}{v(x)u(x)w(x)}dx. $$ Now, suppose that $$ \mathfrak{B}=\operatorname{span}\{B_{0}(x), B_{1}(x),\dots, B_{N}(x)\}, $$ where $\mathfrak{B}$ is a finite-dimensional subspace of $L^{2}_{w}(\Gamma)$, $\operatorname{dim}\mathfrak{B} = N+1$, so $\mathfrak{B}$ is a closed subspace of $L^{2}(\Gamma)$. Therefore, $\mathfrak{B}$ is a complete subspace of $L^{2}(\Gamma)$. Assume that $f(x)$ is an arbitrary element in $L^{2}(\Gamma)$. Thus $f$ has a unique best approximation in $\mathfrak{B}$ subspace, say $\hat{b}(x)\in \mathfrak{B}$; that is, \begin{equation} \exists~ \hat{b}(x)\in\mathfrak{B}, \quad \forall b(x)\in \mathfrak{B},\quad \| f(x)-\hat{b}(x)\| \leq \| f(x)-b(x)\|. \end{equation} Notice that we can write $b(x)$ vector as a combination of the basis vectors of $\mathfrak{B}$ subspace. We know function of $f(x)$ can be expanded by $N+1$ terms of Bessel polynomials as: \[ f(x)=f_{N}(x)+R(x); \] that is, \begin{equation}\label{eqqq07} f_{N}(x)=\sum^{N}_{n=0}{a_{n}B_{n}(x)}=A^{T}B(x), \end{equation} where $B(x)=[B_{0}(x), B_{1}(x),\dots, B_{N}(x)]^{T}$ and $R(x)\in\mathfrak{B}^{\perp}$ that ${\mathfrak B}^{\perp}$ is the orthogonal complement. So $(f(x)-f_{N}(x))\in\mathfrak{B}^{\perp}$ and $b(x)\in\mathfrak{B}$ are orthogonal which we denote it by \[ (f(x)-f_{N}(x))\perp b, \] thus $f(x)-f_{N}(x)$ vector is orthogonal over all of basis vectors of $\mathfrak{B}$ subspace as: \[ \langle f(x)-f_{N}(x),B_{i}(x)\rangle_{w} =\langle f(x)-A^{T}B(x),B_{i}(x)\rangle_{w}=0,\quad i=0, 1,\dots, N, \] hence \[ \langle f(x)-A^{T}B(x),B^{T}(x)\rangle_{w}=0, \] therefore $A$ can be obtained by \begin{gather*} \langle f(x),B^{T}(x)\rangle_{w}=\langle A^{T}B(x),B^{T}(x)\rangle_{w},\\ A^{T}=\langle f(x),B^{T}(x)\rangle_{w}\langle B(x),B^{T}(x)\rangle_{w}^{-1}, \quad n=0, 1,\dots, N. \end{gather*} \subsection{Definition of FRB} Some researchers have proposed the series expansions $\sum_{i=0}^{N}{c_{i}x^{i\alpha}},~(\alpha>0)$ to solve the fractional differential equations, for instance, Bhrawy et al constructed shifted fractional-order Jacobi orthogonal functions to solve the nonlinear initial value problem of fractional order $\alpha$ and a class of time-fractional partial differential equations with variable coefficient \cite{reffff14}. Authors \cite{reffff16,reffff17} have proposed fractional-order Legendre functions to solve fractional-order differential equations and the time-fractional convection-diffusion equation. Alshbool et al. have utilized operational matrices of new fractional Bernstein functions for approximating solutions to fractional differential equations \cite{reffff18}. Parand and Delkhosh have introduced the fractional order of the Chebyshev functions for solving Volterra's population growth model of arbitrary order \cite{reffff188}. Baker \cite{refd07} proved that the answer to Thomas-Fermi equation is as fractional forms, for this reason, we have applied new FRB to solve the Thomas-Fermi equation in the semi-infinite interval, $\{ FB_{n}\}$: \[ FB_{n}^{\alpha}(x,L)=B_{n}(\frac {x^{\alpha}}{x^{\alpha} + L}),\quad n=0, 1,\dots,N \] or \begin{equation}\label{eqqq09} FB_{n}^{\alpha}(x,L)=\sum_{r=0}^{[\frac{N-n}{2}]} \frac{(-1)^r}{r!(n+r)!}(\frac{x^{\alpha}}{2(x^{\alpha} + L)})^{2r+n}, \quad n=0, 1,\dots,N \end{equation} where $\alpha>0$, $x\in [0,\infty)$, $B_{n}(x)$ is Bessel polynomials of order $n$, and the constant parameter $L>0$ is a scaling/stretching factor. Let $\Lambda =\{x:0\leq x < \infty \}$ and $L^{2}_{w}(\Lambda )=\{z :\Lambda \to \mathbb{R}| z$ is measurable and $\| z \|_{w} < \infty \}$, where $$ \| z \|_{w}=\Big(\int^{\infty}_{0}|z(x)|^{2}w(x,L)dx\Big)^{1/2}, $$ with $w(x,L)=\frac{\alpha x^{\alpha-1}L}{(x^{\alpha}+L)^{2}}$, is the norm induced by the inner product of the space $L^{2}_{w}(\Lambda)$ as follows: $$ \langle z,g\rangle_{w}=\int^{\infty}_{0}{z(x)g(x)w(x,L)}dx. $$ Now, suppose that $$ \mathfrak{FB}=\operatorname{span}\{FB^{\alpha}_{0}(x,L), FB^{\alpha}_{1}(x,L), \dots, FB^{\alpha}_{N}(x,L)\}, $$ Let $y(x)\in L^{2}(\Gamma)$ be a function defined over interval $[0,\infty)$ can be expanded by $N+1$ terms of FRB as: \begin{equation}\label{eqqq10} y_{N}(x)=\sum^{N}_{n=0}{a_{n}FB^{\alpha}_{n}(x,L)}=A^{T}FB(x,L), \end{equation} where $FB(x,L)=[FB^{\alpha}_{0}(x,L), FB^{\alpha}_{1}(x,L),\dots, FB^{\alpha}_{N}(x,L)]^T$. Hence \begin{equation}\label{eqqq11} \langle y(x)-A^{T} FB(x,L), FB^{T}(x,L)\rangle_{w}=0, \end{equation} therefore $A$ can be obtained by \begin{gather*} \langle y(x), FB^{T}(x,L)\rangle_{w} =\langle A^{T} FB(x,L), FB^{T}(x,L)\rangle_{w},\\ A^{T}=\langle y(x), FB^{T}(x,L)\rangle_{w}\langle FB(x,L), FB^{T}(x,L)\rangle_{w}^{-1},~n=0, 1,\dots, N. \end{gather*} \section{Quasilinearization method (QLM)}\label{sec3} The QLM is a generalization of the Newton-Raphson method \cite{refff01,refff02} to solve the nonlinear differential equation as a limit of approximating the nonlinear terms by an iterative sequence of linear expressions. Bellman and Kalaba have introduced the QLM method about fifty years ago \cite{refff03,refff04}. The QLM techniques are based on the linearization of the high order ordinary/partial differential equation and require the solution of a linear ordinary differential equation at each iteration. Mandelzweig and Tabakin \cite{refff05} have determined general conditions for the quadratic, monotonic and uniform convergence of the QLM method to solve both initial and boundary value problems in nonlinear ordinary $n$-th order differential equations in $N$-dimensional space. Recently, the QLM method has been successfully applied by researchers to solve the various types of fractional differential equations and some ordinary nonlinear equation \cite{refff06,refff07,refff08,refff09}. We have considered second-order nonlinear ordinary differential equations in one variable on the interval $[0, \infty)$ as follows: \begin{equation}\label{eqqq12} \frac{d^2u}{dx^2}=F(u'(x),u(x),x), \end{equation} with the boundary conditions: $u(0)=A$, $u(\infty)=B$, where $A$ and $B$ are real constants and $F$ is nonlinear functions. By using the QLM for solving \eqref{eqqq11} determines the $(r+1)$-th iterative approximation $u_{r+1}(t)$ as a solution of the linear differential equation: \begin{equation}\label{eqqq12b} \frac{d^2u_{r+1}}{dx^2}=F(u'_r,u_r,x)+(u_{r+1}-u_r)F_u(u'_r,u_r,x)+(u'_{r+1}-u'_r)F_{u'}(u'_r,u_r,x), \end{equation} with the boundary conditions \begin{equation}\label{eqqq13} u_{r+1}(0)=A,\quad u_{r+1}(\infty)=B, \end{equation} where $r=0,1,2,\dots$ and the functions $F_u = \partial F/ \partial u$ and $F_{u'} = \partial F / \partial u'$ are functional derivatives of functional $F(u'_r,u_r,x)$. \section{Solution of Thomas-Fermi equation by FRBC-QLM}\label{sec4} By utilizing QLM technique on \eqref{eqqq01}, we have \begin{equation}\label{eqqq14} \frac{d^{2}y_{r+1}(x)}{dx^2} - \frac{3}{2\sqrt{x}}(y_{r}(x))^{1/2}y_{r+1}(x) =-\frac{1}{2\sqrt{x}}(y_{r}(x))^{3/2}, \end{equation} with the boundary conditions: \begin{equation}\label{eqqq15} y_{r+1}(0)=1,\quad y_{r+1}(\infty)=0, \end{equation} where $r = 0, 1, 2,\dots $. For rapid convergence is actually sufficient that the initial guess be sufficiently best to ensure the smallness of just one of the quantity $q_{r} = k||y_{r+1} - y_{r}||$, where $k$ is a constant independent of $r$. Usually, it is advantageous that $y_{0}(t)$ would satisfy at least one of the boundary conditions \eqref{eqqq15} \cite{refff06}, thus set $y_{0}(x)=1$ for the initial guess of Thomas-Fermi equation. In this paper have been considered two terms $\frac{1}{x^2+1}$ and $\frac{x}{x^2+1}$ to satisfy boundary conditions \eqref{eqqq15}. Thus we can approximate $y_{r+1}(x)$ by $N+1$ basis of FRB as: \begin{equation}\label{eqqq16} y_{r+1}(x) \thickapprox y_{N,r+1}(x) =\frac{1}{x^2+1}+\frac{x}{x^2+1}\sum^{N}_{n=0}{\hat{c_{i}}FB^{\alpha}_{n}(x,L)}. \end{equation} where $\alpha>0$ and $r = 0, 1, 2,\dots $. In all of the spectral methods, the purpose is to find $\hat{c_{i}}$ coefficients. To apply the collocation method, we constructed the residual function for $(r+1)$-th iteration in QLM method by substituting $y_{r+1}(x)$ by $y_{N,r+1}(x)$ into \eqref{eqqq14} as follows: \begin{equation}\label{eqqq17} Res_{r+1}(x)=\frac{d^{2}y_{N,r+1}}{dx^{2}} -\frac{3}{2\sqrt{x}}(y_{r}(x))^{1/2}y_{N,r+1}(x) +\frac{1}{2\sqrt{x}}(y_{r}(x))^{3/2}. \end{equation} A method for forcing the residual function \eqref{eqqq17} to zero can be defined as collocation algorithm. There is no limitation to choose the point in the collocation method. The $N+1$ collocation points which are roots of rational Chebyshev functions on interval $[0,\infty$) (i.e. $x_{i}=(1-cos(\frac{(2i-1)\Pi}{2N+2}))/(1+cos(\frac{(2i-1)\Pi}{2N+2})), i=1, 2,\dots, N+1$ \cite{reff05}) have been substituted $Res_{r+1}(x)$, therefore: \begin{equation}\label{eqqq18} \operatorname{Res}_{r+1}(x_{i})=0,\quad i=0, 1, \dots, N+1. \end{equation} A linear system of equations has been obtained, all of these equations can be solved by Newton method for the unknown coefficients. We have also done all of the computations by Maple 2015 on PC with CPU Core i5, Windows 7 64bit, and 8GB of RAM. Now we can employ the FRBC-QLM iterative algorithm to solve Thomas-Fermi equation as follows: \noindent BEGIN \begin{itemize} \item[] Input variable of $I$ that is the number of iterations of QLM method. \item[] Input variable of $N$ that is the number of basic of the FRB. \item[] Set $y_{N,0}(x)=1$. \end{itemize} \quad For $r=0$ to $I$ do \begin{itemize} \item[] Construct the series \eqref{eqqq16} for approximating $y_{r+1}(x)$ as $y_{N,r+1}(x)$. \item[] Construct the linear differential equation \eqref{eqqq17} by using QLM method on \eqref{eqqq01}. \item[] Substitute $y_{N,r+1}(x)$ into the equation \eqref{eqqq17} and create residual function $Res_{r+1}(x)$.\\ Now we have $N +1$ unknown $\{ \hat{c_{i}}\}_{0}^{N}$. To obtain these unknown coefficients, we need $N +1$ equations. \item[] Choose the roots of order $N+1$ of Rational Chebyshev functions as $N+1$ collocation points: $\{ x_{i}\}_{0}^{N}$. \item[] Substitute collocation points $\{ x_{i}\}_{0}^{N}$ into the $Res_{r+1}(x)$ and create the $N+1$ equations. \item[] Solve the $N+1$ linear equations with $N+1$ unknown coefficients, for calculating $y_{N,r+1}(x)$. \end{itemize} \quad End For\\ END \section{Numerical Results}\label{sec5} The initial slope $y'(0)$ is difficult to compute by any means and plays an important role in determining many physical properties of the Thomas-Fermi atom. It determines the energy of a neutral atom in the Thomas-Fermi approximation. Zaitsev et al \cite{refh6} have shown that the methods of Runge-Kutta and Adams-Bashforth can apply to solve the Thomas-Fermi equation in the semi-infinite interval, although their methods are ill-condition and have not high accuracy for more scheme. Exact solution for Thomas-Fermi differential equation, which is defined on the semi-infinite interval and has a singularity at $x=0$ and its boundary condition occurs at infinity, is not available, so approximating this solution is very important. Table 1 shows a list of the number of calculations $y'(0)$ of the Thomas-Fermi potential. As can be seen, some researchers have achieved good results and accuracy. The last three rows show best approximations of $y'(0)$ for various value of $N$ and a fixed value of $L=1$ by the present method which shows that the present solution is highly accurate. Tables 2 and 3 show values obtained of $y(x)$ and $y'(x)$ by the present method respectively, for different values of $N$ and the 45-th iteration. Obviously, Table 4 and 5 present some numerical example to illustrate the accuracy and convergence of our suggested method by increasing the number of points and iterations. It should be mentioned that all calculations are done by software Maple for various values $N$ and iterations. Figure \ref{fig1} shows the resulting graph of Thomas-Fermi equation obtained by the present method for $N=200$ and iteration 45 which tends to zero as $x$ increases by boundary condition $y(\infty) = 0$, and graphs of residual error of the problem with $N=50,~100,~150,~200$, and the 45-th iteration, note that the residual error decreases with the increase of the collocation points. Comparing the computed results by this method with the others shows that this method provides more accurate and numerically stable solutions than those obtained by other methods. \begin{figure}[ht] \begin{center} \includegraphics[width=6cm]{fig1a} % residual.eps \includegraphics[width=6cm]{fig1b} \\ Graphs of residual error\qquad \hfil Graph of $y(x)$ \end{center} \caption{Graphs of residual error with $N=50,100,150,200$, and iteration 45, and Thomas-Fermi graph obtained by present method.} \label{fig1} \end{figure} \begin{table}[ht] \caption{Comparison of the obtained values of $y'(0)$ by researchers, inaccurate digits are in bold face.} \footnotesize \begin{center} \begin{tabular}{ll} \hline Author/Authors & Obtained value of $y'(0)$ \\ \hline Fermi (1928) \cite{refd05} & -1.58 \\ Baker (1930) \cite{refd07} & -1.588\textbf{558} \\ Bush and Caldwell (1931) \cite{reff11} & -1.58\textbf{9} \\ Miranda (1934) \cite{reff12} & -1.5880\textbf{464} \\ Slater and Krutter (1935) \cite{refd03} & -1.5880\textbf{8} \\ Feynman et al (1949) \cite{refd04} & -1.588\textbf{75} \\ Kobayashi et al. (1955) \cite{refd31} & -1.58807\textbf{0972} \\ Mason (1964) \cite{refd09} & -1.5880710 \\ Laurenzi (1990) \cite{refd06} & -1.588\textbf{588} \\ MacLeod (1992) \cite{reff03} & -1.5880710226 \\ Wazwaz (1999) \cite{refh3} & -1.58807\textbf{6779} \\ Epele et al (1999) \cite{refm18} & -1.588\textbf{1} \\ Esposito (2002) \cite{refd14} & -1.588 \\ Liao (2003) \cite{refd16} & -1.58\textbf{712} \\ Khan and Xu (2007) \cite{refm19} & -1.58\textbf{6494973} \\ El-Nahhas (2008) \cite{refd18} & -1.5\textbf{5167} \\ Yao (2008) \cite{refm12} & -1.5880\textbf{04950} \\ Fernandez (2008) \cite{refm17} & -1.588071022611375313 \\ Parand and Shahini (2009) \cite{reff05} & -1.58807\textbf{02966} \\ Marinca and Herianu (2011) \cite{refh4} & -1.5880\textbf{659888} \\ Oulne (2011) \cite{refh7} & -1.5880710\textbf{34} \\ Abbasbandy and Bervillier (2011) \cite{reff10} & -1.588071022611375312718\textbf{9} \\ Zhu et al. (2012) \cite{refh5} & -1.58\textbf{794357} \\ Turkylmazoglu (2012) \cite{refm11} & -1.5880\textbf{1} \\ Zhao et al (2012) \cite{refm13} & -1.5880710226 \\ Parand et al (2013) \cite{reff09} & -1.58807\textbf{0339} \\ Boyd (2013) (with m=600) \cite{reff01} & -1.5880710226113753127186845 \\ % Tavassoli Kajan et al. (2013) \cite{reff04} & -1.58807102261137\textbf{4} \\ Amore et al (2014) \cite{refm16} & -1.588071022611375312718684508 \\ Marinca and Ene (2014) \cite{refd25} & -1.588071\textbf{9992} \\ Bayatbabolghani and Parand(2014)\cite{reff07} & -1.588071 \\ Kilicman et al (2014) \cite{reff02} & -1.588071\textbf{347} \\ Liu and Zhu (2015) \cite{reff08} & -1.58807\textbf{2} \\ This article [N=100]&-1.5880710226113753127 \\ This article [N=150]& -1.5880710226113753127186845\\ This article [N=200]& -1.588071022611375312718684509423 \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[ht] \caption{Values of $y(x)$ for various values of $x$ with iteration 45 and N=200 } \footnotesize \begin{center} \begin{tabular}{lc|lc} \hline $x$ & $y(x)$ & $x$ & $y(x)$ \\ \hline 0.25 & 0.755201465313331276073659062048 & 30 & 2.255836616202855884224e-3 \\ 0.50 & 0.606986383355979909494446070174 & 40 & 1.113635638833368812571e-3 \\ 0.75 & 0.502346846412368627446521794036 & 50 & 6.322547829849047267797e-4 \\ 1.00 & 0.424008052080705600224612007418 & 60 & 3.939113666854170020266e-4 \\ 1.25 & 0.363201414459514114681451617277 & 70 & 2.622652998120112937417e-4 \\ 1.50 & 0.314777463700458172973580939810 & 80 & 1.835457597407102401554e-4 \\ 1.75 & 0.275451327996091785680113852782 & 90 & 1.335458289537346235905e-4 \\ 2.00 & 0.243008507161119555299806749733 & 100 & 1.002425681394073316855e-4 \\ 2.25 & 0.215894626576130144431137496637 & 200 & 1.450180349694576468040e-5 \\ 2.50 & 0.192984123458000701287136925252 & 300 & 4.548571953616680184257e-6 \\ 2.75 & 0.173441292490063451179594691770 & 400 & 1.979732628112504742575e-6\\ 3.00 & 0.156632673216495841339813440477 & 500 & 1.034077168199939706035e-6 \\ 3.25 & 0.142069642692650781317847467819 & 600 & 6.068687696675251337105e-7 \\ 3.5 & 0.129369596993799111381550504704 & 700 & 3.861765157037986162218e-7 \\ 3.75 & 0.118229001616846686506107664622 & 800 & 2.608137304998336353315e-7 \\ 4.00 & 0.108404256918907711089847680321 & 900 & 1.843724151350651789179e-7 \\ 4.25 & 0.099697845864740046922595440377 & 1000 & 1.351274773541058315394e-7 \\ 4.50 & 0.091948133826563845114632113892 & 2000 & 1.733984751613850910820e-8 \\ 4.75 & 0.085021743728059499429131697623 & 3000 & 5.189408334543857341875e-9 \\ 5.00 & 0.078807779251369904256091892542 & 4000 & 2.201209082423027362721e-9 \\ 6.00 & 0.059422949250422580797949567059 & 5000 & 1.130926706419984771574e-9 \\ 7.00 & 0.046097818604498589876456260102 & 6000 & 6.56056637887703224451e-10 \\ 8.00 & 0.036587255264676802392315804375 & 7000 & 4.13886522087042381058e-10 \\ 9.00 & 0.029590935270546873724362041806 & 8000 & 2.77658195353364145555e-10 \\ 10.00 & 0.024314292988680864190110388176 & 9000 & 1.95225879692802747067e-10 \\ 20.00 & 0.005784941191566940442010571504 &10000& 1.42450044462688615523e-10 \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[ht] \caption{Values of $y'(x)$ for various values of $x$ with iteration 45 and $N=200$} \footnotesize \begin{center} \begin{tabular}{lc|lc} \hline $x$ & $y'(x)$ & $x$ & $y'(x)$ \\ \hline 0.25 &- 0.722306984910234919519668083864 & 30 & -1.806700064769926350e-4 \\ 0.50 & -0.489411612574538088647005847557 & 40 & -6.966802854032586631e-5 \\ 0.75 & -0.358306880167513621987250767311 & 50 & -3.249890204825881462e-5 \\ 1.00 & -0.273989051593306251989464686519 & 60 & -1.719770008309986259e-5 \\ 1.25 & -0.215794130300733601274300529492 & 70 & -9.956533393052361268e-6 \\ 1.50 & -0.173738799013945185681936465228 & 80 & -6.166195528764075475e-6 \\ 1.75 & -0.142320937196893658885960452965 & 90 & -4.024473703766734693e-6 \\ 2.00 & -0.118243191625487620571255867534 & 100 & -2.739351068678330086e-6 \\ 2.25 & -0.099409321201447030009355248089 & 200 & -2.057532316475268926e-7 \\ 2.50 & -0.084426186798809043812545918214 & 300 & -4.365949618530290454e-8 \\ 2.75 & -0.072335044846097235621822505698 & 400 & -1.436682305996181021e-9 \\ 3.00 & -0.062457130854120976228704899999 & 500 & -6.034363442475256759e-9 \\ 3.25 & -0.054300422911798016711579461695 & 600 & -2.961822515102276227e-9 \\ 3.5 & -0.047501046582295208950097689053 & 700 & -1.619832187577198029e-9 \\ 3.75 & -0.041785207716826396797995883443 & 800 & -9.59243855994648160e-10 \\ 4.00 & -0.036943757824123486354813738987 & 900 & -6.03766177055436240e-10 \\ 4.25 & -0.032814785443993540385309640573 & 1000 & -3.98801070822799359e-10 \\ 4.50 & -0.029271448448803843379269452384 & 2000 & -2.57608536991971054e-11 \\ 4.75 & -0.026213311168397937715703460851 & 3000 & -5.15300117640232003e-12 \\ 5.00 & -0.023560074954700512881180449080 & 4000 & -1.64161860704260568e-12 \\ 6.00 & -0.015867549533407079812737615662 & 5000 & -6.75339712187831946e-13 \\ 7.00 & -0.011142531814867088405578460800 & 6000 & -3.26676807336313998e-13 \\ 8.00 & -0.008088602969645474322126751847 & 7000 & -1.76730816571737264e-13 \\ 9.00 & -0.006033074714457392439143608348 & 8000 & -1.03777992089730866e-13 \\ 10.00 & -0.004602881871269254502543511851 & 9000 & -6.48790244833915206e-14 \\ 20.00 & -0.000647254332777692033047085589 &10000& -4.26161649603483992e-14 \\ \hline \end{tabular} \end{center} \end{table} %%%%%%%%%%% \begin{table}[ht] \caption{Numerical solution $y(x)$ with various values of $x$, $N$ and iterations.} \footnotesize \begin{center} \begin{tabular}{| m{0.8em} | m{0.9em} | m{8.3em} m{16.5em} |} %m{15.2em} | } \hline $N$ & $x$ &\quad 15th iteration & \quad 45th iteration\\ \hline % n=50 50 &10 100 200 300 400 500 & 0.024314292988800 0.000100242576939 0.000014501554262 0.000004548556111 0.000001982372128 0.000001042194934 & 0.024314292988680865622793862360 0.000100242568139361977721849610 0.000014501803498894699818187404 0.000004548571957423423877537837 0.000001979732627197689616236145 0.000001034077156421974180242224\\ \hline \hline % n=100 100 & 10 100 200 300 400 500 & 0.024314292988681 0.000100242573102 0.000014501912109 0.000004549264630 0.000001982297142 0.000001040782592& 0.024314292988680864190110392609 0.000100242568139407331736524506 0.000014501803496945768054100507 0.000004548571953616663090340995 0.000001979732628112377070251772 0.000001034077168200025910092884 \\ \hline \hline % n=150 150& 10 100 200 300 400 500 & 0.024314292988682 0.000100242577603 0.000014502011929 0.000004549891945 0.000001984616067 0.000001046957055 & 0.024314292988680864190110388161 0.000100242568139407331685518495 0.000014501803496945764680397208 0.000004548571953616680172305663 0.000001979732628112504797936584 0.000001034077168199940285585132\\ \hline \hline %n=200 200&10 100 200 300 400 500 & 0.024314292988682 0.000100242578682 0.000014502034840 0.000004550038335 0.000001985153638 0.000001048360910 & 0.024314292988680864190110388176 0.000100242568139407331685585932 0.000014501803496945764680403612 0.000004548571953616680184257373 0.000001979732628112504742575794 0.000001034077168199939706035334 \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[ht] \caption{Numerical solution $y'(x)$ with various values of $x$, $N$ and iterations.} \footnotesize \begin{center} \begin{tabular}{| m{1.5em} | m{1.5em} | m{8.2em} m{16.5em} | } \hline $~N$ & $x$ &\quad 15th iteration &\quad 45th iteration\\ \hline % n=50 50 &0 10 100 200 300 400 500 & -1.58798412034597 -0.00460288187129 -0.00000273935089 -0.00000020575729 -0.00000004364704 -0.00000001432569 -0.00000000596694& -1.588071022461220318498896590154 -0.004602881871269254843810395419 -0.000002739351068678834237868374 -0.000000205753231612157754759952 -0.000000043659496195672746505803 -0.000000014366823142132220490812 -0.000000006034363572464639504768\\ \hline \hline % n=100 100 &0 10 100 200 300 400 500 & -1.58806849943926 -0.00460288187126 -0.00000273935085 -0.00000020575078 -0.00000004364886 -0.00000001433829 -0.00000000597884& -1.588071022611375312724621425500 -0.004602881871269254502543510545 -0.000002739351068678330084039598 -0.000000205753231647526755494256 -0.000000043659496185303784906289 -0.000000014366823059962315807005 -0.000000006034363442469561014542\\ \hline \hline % n=150 150&0 10 100 200 300 400 500 & -1.58807102261138 -0.00460288187126 -0.00000273935065 -0.00000020574848 -0.00000004365807 -0.00000001431185 -0.00000000592789 & -1.588071022611375312718684517975 -0.004602881871269254502543511856 -0.000002739351068678330086729788 -0.000000205753231647526892514803 -0.000000043659496185302905165385 -0.000000014366823059961806903847 -0.000000006034363442475252124155 \\ \hline \hline %n=200 200 &0 10 100 200 300 400 500 & -1.58807102261137 -0.00460288187126 -0.00000273935060 -0.00000020574799 -0.00000004363715 -0.00000001430608 -0.00000000591587 & -1.588071022611375312718684509423 -0.004602881871269254502543511851 -0.000002739351068678330086744132 -0.000000205753231647526892605535 -0.000000043659496185302904545958 -0.000000014366823059961810213641 -0.000000006034363442475256759610 \\ \hline \end{tabular} \end{center} \end{table} \section{Conclusion}\label{sec6} The fundamental goal of this paper has been to construct an approximation to the solution of nonlinear Thomas-Fermi equation in a semi-infinite interval which has a singularity at $x = 0$ and its boundary condition occurred in infinity. In the above discussion, we applied a new method to solve the Thomas-Fermi equation that is nonlinear ordinary differential equation on a semi-infinite interval. By using an analytical method for solving Thomas-Fermi equation has proved that the answer to this problem is as fractional forms \cite{refd07}. So for the first time, we solved the problem based on the new fractional order of rational Bessel functions without any domain truncation or transformation of the domain of the problem to a finite domain. In this work, first, by utilizing QLM over Thomas-Fermi equation a sequence of linear differential equations is obtained. Second, at each iteration the linear differential equation is solved by novel FRBC method. We obtained accurately to 30 decimal places for initial slope, $y'(0) = -1.588071022611375312718684509423$, only by using 200 collocation points and successfully have been applied to find the most accurate values of $y(x)$ and $y'(x)$. A known open problem in spectral methods is finding the optimal value for $L$ \cite{reff001}, but in this paper, for simplicity, we set $L = 1$. The numerical results of solving this problem show that this method is higher accurate than obtained results of other famous methods. Finally, the comparison results have shown that the present method is an acceptable approach and good candidate to solve this type of problems that occur in the semi-infinite interval and the nonlinear singular two point boundary value problems effectively. \begin{thebibliography}{00} \bibitem{reff10} S. Abbasbandy, C. Bervillier; \emph{Analytic continuation of Taylor series and the boundary value problems of some nonlinear ordinary differential equations}, Appl. Math. Comput., 218 (2011), 2178-2199. \bibitem{reffff17} S. Abbasbandya, S. Kazemb, M. S. Alhuthali, H. H. Alsulami; \emph{Application of the operational matrix of fractional order Legendre functions for solving the time-fractional convection-diffusion equation}, Appl. Math. Model., 266 (2015), 31-40. \bibitem{refh1} G. Adomian; \emph{Solution of the Thomas-Fermi Equation}, Appl. Math. Lett., 11 (1998), 131-133. \bibitem{reffff18} M. H. T. Alshbool, A. S. Bataineh, I. Hashim, O. R. Isik; \emph{Solution of fractional-order differential equations based on the operational matrices of new fractional Bernstein functions}, J. King Saud Uni. Sci. (2015) doi:10.1016/j.jksus.2015.11.004. \bibitem{refm16} P. Amore, J. P. Boyd, F. M. Fernandez; \emph{Accurate calculation of the solutions to the Thomas-Fermi equations}, Appl. Math. Comput., 232 (2014), 929-943. \bibitem{refm16.5} P. Amore, J. P. Boyd, F. M. Fernandez; \emph{Accurate calculation of the solutions to the Thomas-Fermi equations}, arXiv:1205.1704v2, 2014. \bibitem{refd07} E. B. Baker; \emph{The application of the Fermi-Thomas statistical model to the calculation of potential distribution in positive ions}, Quart. Appl. Math., 36 (1930), 630-647. \bibitem{reff07} F. Bayatbabolghani, K. Parand; \emph{Using Hermite Function for Solving Thomas-Fermi Equation}, Int. J. Math., Comput., Phys, Electr. Comput. Eng., 8(1) (2014), 123-126. \bibitem{reffff03} W. W. Bell; \emph{Special functions for scientists and engineers, D. Van Nostrand Company}, CEf Canada, 1967. \bibitem{refff04} R. E. Bellman, R. E. Kalaba; \emph{Quasilinearization and Nonlinear Boundary-Value Problems, Elsevier Publishing Company}, New York, 1965. \bibitem{reffff14} A. Bhrawy, M. A. Zakyb; \emph{A fractional-order Jacobi Tau method for a class of time-fractional PDEs with variable coefficients}, Math. Method Appl. Sci., 39 (2016), 1765-1779. \bibitem{reff001} J. P. Boyd; \emph{Chebyshev and Fourier spectral Methods}, Second Edition 2000. \bibitem{reff01} J. P. Boyd; \emph{Rational Chebyshev series for the Thomas-Fermi function: Endpoint singularities and spectral methods}, J. Comput. Appl. Math., 244 (2013), 90-101. \bibitem{refdd2} W. Bu, Y. Ting, Y. Wu, J. Yang; \emph{Finite difference/finite element method for two-dimensional space and time fractional blochtorrey equations}, J. Comput. Phys., 293 (2015), 264-279. \bibitem{reff11} V. Bush, S. H. Caldwell; \emph{Thomas-Fermi equation solution by the differential analyzer}, Phys. Rev., 38 (1931), 1898-1902. \bibitem{refd02} S. Chandrasekhar; \emph{Introduction to the Study of Stellar Structure}, Dover, New York, 1967. \bibitem{refdd3} H. J. Choi, J. R. Kweon; \emph{A finite element method for singular solutions of the Navier-Stokes equations on a non-convex polygon}, J. Comput. Appl. Math., 292 (2016) 342-362. \bibitem{reffff01} M. P. Coleman; \emph{An introduction to partial differential equations with MATLAB}, Second Edition 2013. \bibitem{refff01} S. D. Conte, C. de Boor; \emph{Elementary Numerical Analysis}, McGraw-Hill International Editions, 1981. \bibitem{refdd12} O. Coulaud, D. Funaro, O. Kavian; \emph{Laguerre spectral approximation of elliptic problems in exterior domains}, Comput. Methods Appl. Mech. Eng., 80 (1990) 451-458. \bibitem{refd28} C.A. Coulson, N.H. March; \emph{Momenta in Atoms using the Thomas-Fermi Method}, Proc. Phys. Soc., Sect. A, 63(4) (1949), 67-374. \bibitem{refd01} H. T. Davis; \emph{Introduction to Nonlinear Differential and Integral Equations}, Dover, New York, 1962. \bibitem{reffff09} M. Delkhosh; \emph{The conversion a Bessel's equation to a self-adjoint equation and applications}, World Appl. Sci. J., 15 (2011), 1687-1691. \bibitem{refff09} J. V. Devi, F. A. McRae, Z. Drici; \emph{Generalized quasilinearization for fractional differential equations}, Comput. Math. Appl. 59 (2010), 1057-062. \bibitem{refdd16} E. H. Doha, A. H. Bhrawy, R. M. Hafezd, R. A. Gorder; \emph{Jacobi rational-Gauss collocation method for Lane-Emden equations of astrophysical significance}, Nonlinear Anal. Model. Control, 19 (2014), 537-550. \bibitem{refd18} A. El-Nahhas; \emph{Analytic Approximations for Thomas-Fermi Equation}, Acta Phys. Pol. A, 114(4) (2008), 913-918. \bibitem{refm18} L. N. Epele, H. Fanchiotti, C. A. G. Canal, J. A. Ponciano; \emph{Pad\'{e} approximate approach to the Thomas-Fermi problem}, Phys. Rev. A, 60 (1999), 280-283. \bibitem{refd14} S. Esposito; \emph{Majorana solution of the Thomas-Fermi equation}, Am. J. Phys., 70 (2002), 852-856. \bibitem{refd05} E. Fermi; \emph{Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente}, Z. Phys., 48 (1928), 73-79. \bibitem{refm17} F. M. Fernandez; \emph{Rational approximation to the Thomas-Fermi equations}, Appl. Math. Comput., 217 (2011), 6433-6436. \bibitem{refd04} R. P. Feynman, N. Metropolis, E. Teller; \emph{Equations of state of elements based on the generalized Fermi-Thomas theory}, Phys. Rev., 75 (1949), 1561-1573. \bibitem{refdd9} B. Y. Guo; \emph{Error estimation of Hermite spectral method for nonlinear partial differential equations}, Math. Comput., 68 (1999), 1067-1078. \bibitem{Guo1} B. Y. Guo; \emph{Gegenbauer Approximation and Its Applications to Differential Equations on the Whole Line}, J. Math. Anal. Appl., 226 (1998), 180-206. \bibitem{Guo2} B. Y. Guo; \emph{Jacobi Approximations in Certain Hilbert Spaces and Their Applications to Singular Differential Equations}, J. Math. Anal. Appl., 243 (2000), 373-408. \bibitem{refdd19} I. Hashim, M. S. M. Noorani, M. R. S. Hadidi; \emph{Solving the generalized Burgers-Huxley equation using the Adomian decomposition method}, Math. Comput. Model., 43 (2006) 1404-1411. \bibitem{refdd21} J. H. He; \emph{Homotopy perturbation technique}, Comput. Methods in Appl. Mech. Eng., 178 (1999) 257-262. \bibitem{refdd23} J. H. He, X. H. Wu; \emph{Exp-function method for nonlinear wave equations}, Chaos, Soliton. Fract., 30 (2006), 700-708. \bibitem{reffff02} R. L. Herman; \emph{A Course in Mathematical Methods for Physicists}, 2013. \bibitem{refdd17} S. A. Hossayni, J. A. Rad, K. Parand, S. Abbasbandy; \emph{Application of the exact operational matrices for solving the Emden-Fowler equations, arising in? Astrophysics}, Int. J. Ind. Math., 7 (2015) 351-374. \bibitem{refdd15} O. R. Isik, M. Sezer, Z. Guney; \emph{A rational approximation based on Bernstein polynomials for high order initial and boundary values problems}, Appl. Math. Comput., 217 (2011), 9438-9450. \bibitem{reff06} R. Jovanovic, S. Kais, F. H. Alharbi; \emph{spectral Method for Solving the Nonlinear Thomas-Fermi Equation Based on Exponential Functions}, J. Appl. Math., (2014) Article ID 168568, 8 pages. \bibitem{reffff16} S. Kazem, S. Abbasbandy, S. Kumar; \emph{Fractional-order Legendre functions for solving fractional-order differential equations}, Appl. Math. Model., 37 (2013), 5498-5510. \bibitem{refff03} R. Kalaba; \emph{On nonlinear differential equations, the maximum operation and monotone convergence, RAND Corporation}, P-1163, (1957). \bibitem{refm19} H. Khan, H. Xu; \emph{Series solution to the Thomas-Fermi equation}, Physics Letters A, 365 (2007), 111-115. \bibitem{reff02} A. Kilicman, I. Hashimb, M. Tavassoli Kajani, M. Maleki; \emph{On the rational second kind Chebyshev pseudospectral method for the solution of the Thomas-Fermi equation over an infinite interval}, J. Comput. Appl. Math., 257 (2014), 79-85. \bibitem{refd31} S. Kobayashi, T. Matsukuma, S. Nagi, K. Umeda; \emph{Accurate value of the initial slope of the ordinary T-F function, J. Phys. Soc. Japan}, 10 (1955), 759-762. \bibitem{refd06} B. J. Laurenzi; \emph{An analytic solution to the Thomas-Fermi equation}, J. Math. Phys., 10 (1990) 2535-2537. \bibitem{reff08} C. Liu, S. Zhu; \emph{Laguerre pseudospectral approximation to the Thomas-Fermi equation}, J. Comput. Appl. Math., 282 (2015) 251-261. \bibitem{refd16} S. Liao, \emph{An explicit analytic solution to the Thomas-Fermi equation}, Appl. Math. Comput., 144 (2003) 495-506. \bibitem{reff03} A. J. MacLeod; \emph{Chebyshev series solution of the Thomas-Fermi equation}, Comput. Phys. Commun., 67 (1992), 389-391. \bibitem{refff05} V. B. Mandelzweig, F. Tabakin; \emph{Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs}, Comput. Phys. Commun. 141 (2001), 268-281. \bibitem{refh4} V. Marinca, N. Herisanu; \emph{An optimal iteration method with application to the Thomas-Fermi equation}, Cent. Eur. J. Phys., 9 (2011), 891-895. \bibitem{refd25} V. Marinca, R.D. Ene; \emph{Analytical approximate solutions to the Thomas-Fermi equation}, Cent. Eur. J. Phys., 12(7) (2014), 503-510. \bibitem{refd09} J. C. Mason; \emph{Rational approximations to the ordinary Thomas-Fermi function and its derivative}, Proc. Phys. Soc., 84 (1964), 357-359. \bibitem{reff12} C. Miranda; \emph{Teoremi e metodi per lintegrazione numerica della equazione differenziale di Fermi}, Memorie della Reale Accademia dItalia, Classe di scienze fisiche, Mat. Nat., 5 (1934), 285-322. \bibitem{refdd1} B. J. Noye, M. Dehghan; \emph{New explicit finite difference schemes for two-dimensional diffusion subject to specification of mass}, Numer. Meth. Part. Diff. Eq., 15 (1999), 521-534. \bibitem{refh7} M. Oulne; \emph{Variation and series approach to the Thomas-Fermi equation}, Appl. Math. Comput., 218 (2011), 303-307. \bibitem{refdd8} K. Parand, A. R. Rezaei, A. Taghavi; \emph{Numerical approximations for population growth model by rational Chebyshev and Hermite functions collocation approach: a comparison}, Math. Method Appl. Sci., 33 (2010), 2076-2086. \bibitem{reff055} K. Parand, H. Yousefi, M. Delkhosh, A. Ghaderi; \emph{A novel numerical technique to obtain an accurate solution to the Thomas-Fermi equation}, Eur. Phys. J. Plus, 131(7) (2016), 228. \bibitem{reffff07} K. Parand, J. A. Rad, M. Nikarya; \emph{A new numerical algorithm based on the first kind of modified Bessel function to solve population growth in a closed system}, Int. J. Comput. Math., 91 (2014), 1239-1254. \bibitem{refff07} K. Parand, M. Ghasemi, S. Rezazadeh, A. Peiravi, A. Ghorbanpour, A. T. Golpaygani; \emph{Quasilinearization approach for solving Volterra's population model}, Appl. Comput. Math., 9 (2010), 95-103. \bibitem{refdd5} K. Parand, M. Hemami; \emph{Application of Meshfree Method Based on Compactly Sup-ported Radial Basis Function for Solving Unsteady Isothermal Gas Through a Micro-Nano Porous Medium}, Iran. J. Sci. Tech. Trans. A: Science, 1-1, (2015). \bibitem{reff09} K. Parand, M. Dehghan, A. Pirkhedri; \emph{The Sinc-collocation method for solving the Thomas-Fermi equation}, J. Comput. Appl. Math., 237 (2013), 244-252. \bibitem{refdd11} K. Parand, M. Dehghan, A. Taghavi; \emph{Modified generalized Laguerre function Tau method for solving laminar viscous flow The Blasius equation}, Int. J. Numer. Method. H., 20 (2010), 728-743. \bibitem{refdd24} K. Parand, M. Dehghan, F. Baharifard; \emph{Solving a laminar boundary layer equation with the rational Gegenbauer functions}, Appl. Math. Model., 37 (2013), 851-863. \bibitem{refdd6} K. Parand, M. Delkhosh; \emph{Numerical Solution of an Integro-Differential Equation Arising in Oscillating Magnetic Fields}, J. Korean Soc. Indus. Appl. Math., 20 (3), 261-275. \bibitem{reffff188} K. Parand, M. Delkhosh; \emph{Solving Volterra's population growth model of arbitrary order using the generalized fractional order of the Chebyshev functions}, Ricerche Mat., 65(1) (2016), 307-328. \bibitem{refdd10} K. Parand, M. Delkhosh; \emph{Solving the nonlinear Schlomilch's integral equation arising in ionospheric problems}, Afr. Mat., (2016) doi:10.1007/s13370-016-0459-3. \bibitem{reffff08} K. Parand, M. Nikarya; \emph{Solving the Unsteady Isothermal Gas Through a Micro-Nano Porous Medium via Bessel Function Collocation Method}, J. Comput. Theor. Nanosci., 11 (2014), 131-136. \bibitem{reffff04} K. Parand, M. Nikarya, J. A. Rad, F. Baharifard; \emph{A new reliable numerical algorithm based on the first kind of Bessel functions to solve Prandtl-Blasius laminar viscous flow over a semi-infinite flat plate}, Z. Naturforsch. A, 67 (2012), 665-673. \bibitem{reff05} K. Parand , M. Shahini; \emph{Rational Chebyshev pseudospectral approach for solving Thomas-Fermi equation}, Phys. Lett. A, 373 (2009), 210-213. \bibitem{refdd25} K. Parand, A. Taghavi, M. Shahini; \emph{Comparison between rational Chebyshev and modified generalized Laguerre functions pseudospectral methods for solving Lane-Emden and unsteady gas equations}, Acta Phys. Pol. B, 40(6) (2009), 1749-1763. \bibitem{refdd4} K. Parand, S. Abbasbandy, S. Kazem, A .R. Rezaei; \emph{An improved numerical method for a class of astrophysics problems based on radial basis functions}, Phys. Scripta, 83(1) (2011), 015011. \bibitem{refdd5_1} K Parand, S Hashemi; \emph{RBF-DQ Method for Solving Non-linear Differential Equations of Lane-Emden type}, Ain Shams Engin. J., (2016) doi: 10.1016/j.asej.2016.03.010. \bibitem{refdd14} K. Parand, Z. Delafkar, N. Pakniat, A. Pirkhedri, M. Kazemnasab Haji; \emph{Collocation method using Sinc and Rational Legendre functions for solving Volterra's population model}, 16 (2011), 1811-1819. \bibitem{refff06} A. Rezaei, F. Baharifard, K. Parand; \emph{Quasilinearization-Barycentric approach for numerical investigation of the boundary value Fin problem}, Int. J. Comput., Electr., Autom., Control Inf. Eng., 5 (2011), 194-201. \bibitem{refff02} A. Ralston, P. Rabinowitz; \emph{A First Course in Numerical Analysis}, McGraw-Hill Inter-national Editions, 1988. \bibitem{refdd211} H. Saeedi, F. Samimi, \emph{He's homotopy perturbation method for nonlinear ferdholm integro-differential equations of fractional order}, Int. J. Eng. Res. Appl., 2(5) (2012), 52-56. \bibitem{reffff10} N. Sahin, S. Yuzbasi, M. Gulsu; \emph{A collocation approach for solving systems of linear Volterra integral equations with variable coefficients}, Comput. Math. Appl., 62 (2011), 755-769. \bibitem{refdd22} F. Shakeri and M. Dehghan; \emph{Numerical solution of the Klein-Gordon equation via He's variational iteration method}, Nonlinear Dynam., 51 (2008) 89-97. \bibitem{refd03} J. C. Slater, H. M. Krutter; \emph{The Thomas-Fermi method for metals}, Phys. Rev., 47 (1935), 559-568. \bibitem{refdd20} M. Tatari, M. Dehghan, M. Razzaghi; \emph{Application of the Adomian decomposition method for the Fokker-Planck equation}, Math. Comput. Model., 45 (2007), 639-650. \bibitem{refd00} L. H. Thomas; \emph{The calculation of atomic fields, Math. Proc. Cambridge Philos}. Soc., 23 (1927), 542-548. \bibitem{reffff1302} E. Tohidi, H. Nik; \emph{A Bessel collocation method for solving fractional optimal control problems}, Appl. Math. Model., 39 (2015), 455-465. \bibitem{refm11} M. Turkyilmazoglu; \emph{Solution of the Thomas-Fermi equation with a convergent approach}, Commun. Nonlinear. Sci. Numer. Simulat., 17 (2012), 4097-4103. \bibitem{refh3} A. M. Wazwaz, \emph{The modified decomposition method and Pad\'{e} approximates for solving the Thomas-Fermi equation}, Appl. Math. Comput., 105 (1999) 11-19. \bibitem{refff08} A. Yakar; \emph{Initial time difference quasilinearization for Caputo fractional differential equations}, Adv. Diff. Eq., 1 (2012), 1-9. \bibitem{refm12} B. Yao; \emph{A series solution to the Thomas-Fermi equation}, Appl. Math. Comput., 203 (2008), 396-401. \bibitem{reffff13} S. Yuzbasi; \emph{A numerical approach for solving a class of the nonlinear Lane-Emden type equations arising in astrophysics}, Math. Method Appl. Sci., 34 (2011), 2218-2230. \bibitem{reffff11} S. Yuzbasi, M. Sezer; \emph{A numerical method to solve a class of linear integro-differential equations with weakly singular kernel}, Math. Method Appl. Sci., 35 (2012), 621-632. \bibitem{reffff12} S. Yuzbasi, N. Sahin, M. Sezer; \emph{Bessel polynomial solutions of high-order linear Volterra integro-differential equations}, Comput. Math. Appl., 62 (2011), 1940-1956. \bibitem{refh6} N. A. Zaitsev, I. V. Matyushkin, D. V. Shamonov; \emph{Numerical Solution of the Thomas-Fermi Equation for the Centrally Symmetric Atom}, Russ. Microlectron., 33 (2014), 372-378. \bibitem{refh5} S. Zhu, H. Zhu, Q. Wu, Y. Khan; \emph{An adaptive algorithm for the Thomas-Fermi equation}, Numer. Algor., 59 (2012), 359-372. \bibitem{refm13} Y. Zhao, Z. Lin, Z. Liu, S. Liao; \emph{The improved homotopy analysis method for the Thomas-Fermi equation}, Appl. Math. Comput., 218 (2012), 8363-8369. \end{thebibliography} \end{document}