In recent papers the technique for a local and global error estimation and the local-global step size control were presented to solve both ordinary differential equations and semi-explicit index 1 differential-algebraic systems by multistep methods with any reasonable accuracy attained automatically. Now those results are extended to the concept of multistep extrapolation, and the paper demonstrates with numerical examples how such methods work in practice. Especially, we develop an efficient technique to calculate higher derivatives of a numerical solution with Hermite interpolating polynomials. The necessary theory is also provided.