New methods are proposed for the numerical solution of systems of first-order differential equations. On each subinterval of a given mesh of size h, a polynomial of degree l is constructed, its parameters being determined by a multiple collocation technique. The resulting piecewise polynomial approximation is of order
O
(
h
l
+
1
)
O({h^{l + 1}})
at the mesh points and between them. In addition, the jth derivatives of the approximation on each subinterval provide approximations of order
O
(
h
l
+
1
−
j
)
O({h^{l + 1 - j}})
,
j
=
1
,
…
,
l
j = 1, \ldots ,l
. Some of the methods proposed are shown to be A-stable or even strongly A-stable.