Systems of Lane–Emden-type equations describe the modelling of catalytic diffusion reactions, the concentration of carbon substrate and oxygen, the steady state concentration of carbon dioxide, dusty fluid models, and pattern formation. In this article, we explore a collocation method that is based on mth-degree shifted Jacobi polynomials \(\textsf{P}_{L,m}^{(\alpha ,\beta )}(t):= P_{m}^{(\alpha ,\beta )}\left( 2t/L-1\right) \) \((t\in [0,L],L>0;\alpha ,\beta \in (-1,\infty ))\) to obtain numerical solutions of systems of Lane–Emden-type equations. The proposed method assumes the solution in the form of shifted Jacobi polynomial series; thus, the differentiation formula for the shifted Jacobi polynomial is applied upon substituting the assumed series solution into the proposed problem. Collocating at the nodes of the shifted Jacobi–Gauss interpolation in the interval [0, L], we obtain a set of nonlinear algebraic equations, which are subsequently solved for the expansion coefficients using Newton’s iteration method. The convergence and error analysis of the proposed collocation scheme is presented. Two special systems of nonlinear singular initial value problems of Lane–Emden-type (with different values of parameters \(\alpha ,\beta \) ) are presented to demonstrate the proposed method’s reliability, effectiveness, and accuracy. The obtained numerical solutions are compared with the exact solutions and other published results. Our results are in excellent agreement with the solutions under comparison, which is a clear indication that the proposed method is highly reliable, effective, and accurate for nonlinear singular initial value problems.