A class of methods for the numerical solution of systems of ordinary differential equations is given which—for linear systems—gives solutions which conserve the stability property of the differential equation. The methods are of a quadrature type $$y_{i,r} = y_{n,r - 1} + h\sum\limits_{k = 1}^n {a_{ik} f(y_{k,r} ), n = 1,2, \ldots ,n, r = 1,2, \ldots ,} y_{n,0} given$$ wherea ik are quadrature coefficients over the zeros ofP n −P n−1 (v=1) orP n −P n−2 (v=2), whereP n is the Legendre polynomial orthogonal on [0,1] and normalized such thatP m (1)=1. It is shown that $$\left| {y_{n,r} - y(rh) = 0(h^{2n - _v } )} \right|$$ wherey is the solution of $$\frac{{dy}}{{dt}} = f(y), t \mathbin{\lower.3ex\hbox{$\buildrel>\over{\smash{\scriptstyle=}\vphantom{_x}}$}} 0, y(0) given.$$