In this work, we introduce a numerical method for solving nonlinear fractional system of Lane–Emden type equations. The proposed technique is based on Dickson operational matrix of a fractional derivative. First, we deduce the Dickson operational matrix of the fractional derivative using Dickson polynomial, and then, the obtained matrix is unitized to convert the fractional Lane–Emden system with its initial conditions into a system of nonlinear algebraic equations. This system of algebraic equations can be solved numerically via Newton’s iteration method. An error estimate of the proposed method is derived. Numerical examples are provided to demonstrate the validity, applicability, and accuracy of the new technique.