A finite-difference time-domain (FDTD) method for simulating wave propagation in Havriliak-Negami (H-N) dispersive media is presented. In an H-N medium, the time-domain polarization relation involves a fractional differential operator, whose analytical evaluation is troublesome. Moreover, if the fractional differential operator is approximated by a sum of fractional time derivatives, the numerical treatment of the polarization relation results in significant memory storage demands, because fractional derivatives are not local. However, by an appropriate approximation of the Grunwald-Letnikov definition of fractional derivatives, an FDTD scheme with reasonable memory requirements is feasible. In particular, the coefficients that appear in the Grunwald-Letnikov definition are approximated by sums of decaying exponentials. As a result, the FDTD scheme is based on recursive equations, while a limited number of auxiliary vectors need to be stored. The proposed scheme has been applied successfully to the simulation of the excitation of an H-N medium by a wideband Gaussian electromagnetic pulse.