SUMMARYTwo numerical methods for solving two-point boundary-value problems associated with systems of first-order nonlinear ordinary differential equations are described. The first method, which is based on Lobatto quadrature, requires four internal function evaluations for each subinterval. It does not need derivatives and is of order h7, where h is the space chop. The second method, which is similar to the first but is based on Lobatto-Hermite quadrature, makes the additional use of derivatives to achieve O( h9) accuracy. Results of computational experiments comparing these methods with other known methods are given.