The authors consider the following neutral delay-differential equations with multiple delays (NMDDE) \[ y'(t)=Ly(t)+\sum^m_{j=1}[M_j y(t-\tau_j)+N_j y'(t-\tau_j)], \quad t\geq 0, \tag{\(*\)} \] where \(\tau>0, L, M_j\) and \(N_j\) are constant complex-value \(d\times d\) matrices. A sufficient condition for the asymptotic stability of NMDDE system (\(*\)) is given. The stability of Butcher’s (A,B,C)-method for systems of NMDDE are studied. In addition, they present a parallel diagonally-implicit iteration Runge-Kutta (PDIRK) method (NPDIRK) for systems of NMDDE, which is easier to be implemented than fully implicit Runge-Kutta methods. The stability of a special class of NPDIRK methods is investigated by analyzing their stability behaviours of the solutions of (\(*\)). Numerical experiments have not be performed.