An algorithm for the numerical solution of boundary value problems for ordinary differential equations based on the collocation method and representation of the solution as an expansion in Chebyshev polynomials is implemented. It is proposed instead of the traditional approach - merging all conditions (differential and boundary) into one system of linear algebraic equations (SLAE) - to switch to a method for solving the problem in several separate stages. First, spectral coefficients are identified that determine the "general" solution of the original problem. The complexity of reducing the SLAE matrix to a diagonal form (in the case of ODE systems with constant coefficients) at this stage is equivalent to the complexity of multiplying the Chebyshev matrix of coefficients by the vector of the right side of the system. At the second stage, account of the boundary conditions selects a "particular" desired solution, uniquely defining the missing coefficients of the desired expansion. The proposed method can be used to model problems in classical mechanics. © 2023 Tomsk State University. All rights reserved.