The two-stage method for solving ODEs is based on the interpolation of the coefficients of the derivative of the solution. Transforms with a spectral integration matrix gives a part of the solution coefficients. The rest are calculated by the boundary conditions. The complexity of the algorithm is equivalent to the complexity of multiplying the matrix by a vector.