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.