In this paper, a new Legendre orthonormal polynomial–based least–squares method(LSM) is provided for fractional integro-differential equations(FIDEs). We first divide the domain into $n$ cells, and in each cell we apply the least–squares algorithm to obtain an approximating solution. The stability and the optimal convergence order in $W_2^2$–norm error are proven. Numerical examples are studied to verify our theoretical discovery. Comparison with the sextic and eighth $C^{3}-$spline method illustrates that our algorithm can obtain a more accurate approximating solution.