  We propose a numerical method for approximating integro-differential
equations arising in age-of-infection epidemic models. The method is based on a
non-standard finite differences approximation of the integral term appearing in
the equation. The study of convergence properties and the analysis of the
qualitative behavior of the numerical solution show that it preserves all the
basic properties of the continuous model with no restrictive conditions on the
step-length $h$ of integration and that it recovers the continuous dynamic as
$h$ tends to zero.
