A new approach to the numerical solution of systems of first-order ordinary differential equations is given by finding local Galerkin approximations on each subinterval of a given mesh of size
h
h
. One step at a time, a piecewise polynomial, of degree
n
n
and class
C
0
{C^0}
, is constructed, which yields an approximation of order
O
(
h
2
n
)
O({h^{2n}})
at the mesh points and
O
(
h
n
+
1
)
O({h^{n + 1}})
between mesh points. In addition, the
j
j
th derivatives of the approximation on each subinterval have errors of order
O
(
h
n
−
j
+
1
)
,
1
≦
j
≦
n
O({h^{n - j + 1}}),1 \leqq j \leqq n
. The methods are related to collocation schemes and to implicit Runge-Kutta schemes based on Gauss-Legendre quadrature, from which it follows that the Galerkin methods are
A
A
-stable.