Introduction. The paper considers the solution of boundary value problems on an interval for linear ordinary differential equations, in which the coefficients and the right-hand side are continuous functions. The conditions for the orthogonality of the residual equation to the coordinate functions are supplemented by a system of linearly independent boundary conditions. The number of coordinate functions m must exceed the order n of the differential equation.

   Materials and Methods. To numerically solve the boundary value problem, a system of linearly independent coordinate functions is proposed on a symmetric interval [−1,1], where each function has a unit Chebyshev’s norm. A modified Petrov-Galerkin method is applied, incorporating linearly independent boundary conditions from the original problem into the system of linear algebraic equations. An integral quadrature formula with twelfth-order error is used to compute the scalar product of two functions.

   Results. A criterion for the existence and uniqueness of a solution to the boundary value problem is obtained, provided that n linearly independent solutions of the homogeneous differential equation are known. Formulas are derived for the matrix coefficients and the coefficients of the right-hand side in the system of linear algebraic equations for the vector expansion of the solution in terms of the coordinate function system. These formulas are obtained for second- and third-order linear differential equations. The modified Bubnov-Galerkin method is formulated for differential equations of arbitrary order.

   Discussion and Conclusions. The derived formulas for the generalized Bubnov-Galerkin method may be useful for solving boundary value problems involving linear ordinary differential equations. Three boundary value problems with second- and third-order differential equations are numerically solved, with the uniform norm of the residual not exceeding 10–11.