It is known that discontinuities may exist in the solution of neutral delay differential equations even though the function is assumed to be continuous along the interval. This problem occurs when the primary discontinuity in the derivatives solution at the initial point propagates to the subsequent points, which results in a secondary discontinuity. As a result, the solution of the neutral delay may no longer be smooth and lead to a larger number of failure steps. This study proposes a block multistep method to deal with the propagation of derivatives discontinuities in neutral delay. The new invention of the numerical approaches by adapting the block multistep method with the Runge–Kutta Fehlberg variable step strategy is developed. The strategies to approximate both retarded and neutral delays and discontinuity tracking equations are performed to maximize the accuracy of the solution. The error analysis is presented by comparing the numerical results with the existing methods to verify the efficiency of the developed approaches. It is demonstrated that the proposed numerical approaches are able to correct the propagation of discontinuities and provide very smooth solutions with accurate results.
