In this paper, we are concerned with numerical methods for linear neutral systems with multiple delays. For delay-dependently stable neutral systems, we ask what conditions must be imposed on linear multi-step methods in order that the numerical solutions display stability property analogous to that displayed by the exact solutions. Combining with Lagrange interpolation, linear multi-step methods can be applied to the neutral systems. Utilizing the argument principle, a sufficient condition is derived for linear multi-step methods with preserving delay-dependent stability. Numerical examples are given to illustrate the main results.