This paper analyzes the coupled system of nonlinear fractional differential equations involving the caputo fractional derivatives of order $\alpha \in (1,2)$ on the interval (0, T). Our method of analysis is based on the reduction of the given system to an equivalent system of integral equations, then the resulting equation is discretized by using a spectral method based on the Legendre polynomials. We have constructed a Legendre spectral collocation method for the coupled system of nonlinear fractional differential equations. The error bounds under the $L^2-$ and $L^{\mathrm{\infty }}-$norms is also provided, then the theoretical result is validated by a number of numerical tests.