In this paper, we consider the numerical method that is proposed and analyzed in [J. Cao and C. Xu, J. Comput. Phys., 238 (2013), pp. 154–168] for the fractional ordinary differential equations. It is based on the so-called block-by-block approach, which is a common method for the integral equations. We extend the technique to solve the nonlinear system of fractional ordinary differential equations (FODEs) and present a general technique to construct high order schemes for the numerical solution of the nonlinear coupled system of fractional ordinary differential equations (FODEs). By using the present method, we are able to construct a high order schema for nonlinear system of FODEs of the order $α$, $α>0$. The stability and convergence of the schema is rigorously established. Under the smoothness assumption $f,g ∈ C^4[0,T]$, we prove that the numerical solution converges to the exact solution with order $3+α$ for $0<α≤1$ and order $4$ for $α>1$. Some numerical examples are provided to confirm the theoretical claims.