The paper presents a theoretical approach to the construction of extrapolation methods for systems of the kind. {Mathematical expression} where L is a general linear differential operator of order k. For ε=0, the discretization schemes are required to be exact and to contain only solutions in the nullspace of L. For ε≠0, the paper studies the construction of methods that permit quadratic extrapolation. In the special case k=2, a new two-step method is obtained that applies to systems of the type {Mathematical expression} where A is a real, symmetric, positive semi-definite matrix. This algorithm might be of use in regular celestial mechanics-apart from any other possible applications. © 1979 Birkhäuser Verlag.