In this paper we consider the numerical solution of fractional differential equations. In particular, we study a step-by-step procedure, defined over a graded mesh, which is based on a truncated expansion of the vector field along the orthonormal Jacobi polynomial basis. Under mild hypotheses, the proposed procedure is capable of getting spectral accuracy. A few numerical examples are reported to confirm the theoretical findings.