Exact numerical schemes have previously been obtained for some linear retarded delay differential equations and systems. Those schemes were derived from explicit expressions of the exact solutions, and were expressed in the form of perturbed difference systems, involving the values at previous delay intervals. In this work, we propose to directly obtain expressions of the same type for the fundamental solutions of linear delay differential equations, by considering vector equations with vector components corresponding to delay-lagged values at previous intervals. From these expressions for the fundamental solutions, exact numerical schemes for arbitrary initial functions can be proposed, and they may also facilitate obtaining explicit exact solutions. We apply this approach to obtain an exact numerical scheme for the first order linear neutral equation x′(t)−γx′(t−τ)=αx(t)+βx(t−τ), with the general initial condition x(t)=φ(t) for −τ≤t≤0. The resulting expression reduces to those previously published for the corresponding retarded equations when γ=0.

In this work, we obtain an expression, given in terms of some special functions, for the exact numerical solution of the second order delay differential equation with general initial condition for . Based on this exact solution, we propose a family of nonstandard finite difference methods that combine high order accuracy and efficient computational properties. We also show that the numerical solutions obtained with the new schemes are consistent with asymptotic stability properties of the exact solutions.