An FDTD scheme for simulating wave propagation inside Cole-Cole, Havriliak-Negami, or Raicu dispersive media is proposed. The main difficulty in FDTD implementations for such media is the fact that the time-domain polarization relations are fractional integrodifferential equations. To circumvent this difficulty, the dispersive susceptibility of the medium is modeled by Pade approximants. It is proved that under certain conditions the Pade approximant is equal to a weighted sum of Debye terms, while the physical restrictions for positive weights and relaxation times are fulfilled with certainty. Hence, a set of first-order auxiliary differential equations, which approximate the original fractional polarization relation, is derived and directly embedded in the FDTD scheme. In contrast to previous works, the derivation of the Debye expansion of the original dispersive permittivity is carried out in a straightforward manner without involving nonlinear optimization techniques. Comparisons between analytical and FDTD results of the reflection and transmission coefficients as well as of the transfer functions inside the dispersive media illustrate the efficiency of the proposed scheme, over wideband frequency domains