This paper deals with the numerical solution of initial valueproblems for systems of neutral differential equations where т ＞ 0, f and φ denote given vector-valued functions. The numerical stability of a linear multistep method is investigated by analysing the solution of the test equations y’(t) = Ay(t) + By(t -т ) + Cy’(t -т ), where A, B and C denote constant complex N × N-matrices, and т ＞ 0. We investigate the properties of adaptation of the linear multistep method and the characterization of the stability region. It is proved that the linear multistep method is NGP-stable if and only if it is A-stable for ordinary differential equations.