The Lane–Emden type equations are employed in the modeling of several phenomena in the areas of mathematical physics and astrophysics. These equations are categorized as non-linear singular ordinary differential equations on the semi-infinite domain $$[0,\infty )$$ . In this research we introduce the Bessel orthogonal functions as new basis for spectral methods and also, present an efficient numerical algorithm based on them and collocation method for solving these well-known equations. We compare the obtained results with other results to verify the accuracy and efficiency of the presented scheme. To obtain the orthogonal Bessel functions we need their roots. We use the algorithm presented by Glaser et al. (SIAM J Sci Comput 29:1420–1438, 2007) to obtain the $$N$$ roots of Bessel functions.