This paper presents a novel physics-informed machine learning approach, designed to approximate solutions to a specific category of Thomas–Fermi differential equations. To tackle the inherent intricacies of solving Thomas–Fermi equations, we employ the quasi-linearization technique, which transforms the original non-linear problem into a series of linear differential equations. Our approach utilizes collocation least-squares support vector regression, leveraging fractional Chebyshev functions for finite interval simulations and fractional rational Chebyshev functions for semi-infinite intervals, enabling precise solutions for linear differential equations across varied spatial domains. These selections facilitate efficient simulations across both finite and semi-infinite intervals. Key contributions of our approach encompass its versatility, demonstrated through successful approximation of solutions for diverse Thomas–Fermi problem types, including those featuring non-local integral terms, Bohr radius boundary conditions, and isolated neutral atom boundary conditions defined on semi-infinite domains. Furthermore, our method exhibits computational efficiency, surpassing classical collocation methods by solving a sequence of positive definite linear equations or quadratic programming problems. Notably, our approach showcases precision, as evidenced by experiments, including the attainment of the initial slope of the renowned Thomas–Fermi equation with an impressive 39-digit precision.