  We present a novel numerical method, called {\tt Jacobi-predictor-corrector
approach}, for the numerical solution of fractional ordinary differential
equations based on the polynomial interpolation and the Gauss-Lobatto
quadrature w.r.t. the Jacobi-weight function
$\omega(s)=(1-s)^{\alpha-1}(1+s)^0$. This method has the computational cost
O(N) and the convergent order $IN$, where $N$ and $IN$ are, respectively, the
total computational steps and the number of used interpolating points. The
detailed error analysis is performed, and the extensive numerical experiments
confirm the theoretical results and show the robustness of this method.
