This paper deals with the numerical solution of initial value problems for systems of neutral differential equations $$y'(t)=f(t,y(t),y(t- \tau ),y'(t- \tau )), t ＞ 0, $$ $$y(t) = φ(t) \ t＜0,$$ where $\tau＞ 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-\tau) + Cy'(t-\tau),$ where $A, B$ and $C$ denote constant complex $N \times N$-matrices, and $\tau ＞ 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.