  We present a method for the calculation of dynamical correlation functions of
quantum impurity systems out of equilibrium using Wilson's numerical
renormalization group. Our formulation is based on a complete basis set of the
Wilson chain and embeds the recently derived algorithm for equilibrium spectral
functions. Our method fulfills the spectral weight conserving sum-rule exactly
by construction. A local Coulomb repulsion $U>0$ is switched on at $t=0$, and
the asymptotic steady-state spectral functions are obtained for various values
of $U$ as well as magnetic field strength $H$ and temperature $T$. These
benchmark tests show excellent agreement between the time-evolved and the
directly calculated equilibrium NRG spectra for finite $U$. This method could
be used for calculating steady-state non-equilibrium spectral functions at
finite bias through interacting nano-devices.
