We propose a non-perturbative method to simulate heterogeneous electron transfer dynamics in systems described by a Newns-Anderson type of model. The coupling between the molecule and the continuum electronic states at the metal or semiconductor surface is represented using a set of effective modes, by employing an exponential expansion of the bath correlation functions. Depending on the nature of the problems, the nuclear degrees of freedom are either treated explicitly using wave functions and density operators or as dissipative modes using the techniques from the hierarchical equations of the motion method. Numerical examples are also presented for applications in problems including (1) photo-induced charge transfer at the molecule-semiconductor interfaces, (2) heterogeneous electron transfer at the molecule-metal interface, and (3) vibrational relaxation on a metal surface.