Relaxation oscillators pose a challenge for numerical integration, due to the presence of fast and slow time scales. Conventional exponential integrators (such as exponential Euler) allow for numerical stability at large time steps, but do a poor job at capturing the limit cycles governing oscillatory behavior unless the time step size is taken very small. In practice, this results in numerical solutions with the wrong amplitude and/or frequency of oscillation. We present a new family of methods that can maintain stability and preserve limit cycles for much larger time step sizes, with no increase in computational effort over exponential Euler. This is illustrated for the Van der Pol oscillator, as well as for the Hodgkin-Huxley model of neuronal dynamics (whose oscillations correspond to neuronal spiking). Although these systems are dissipative rather than conservative, Hamiltonian and symplectic structures play an interesting role.