We have developed a new algorithm for 3D time-domain electromagnetic (EM) modeling, taking full account of induced polarization (IP) and the coupling between EM and IP effects. The algorithm can be used to model grounded source IP surveys that indicate EM induction effects and airborne time-domain EM surveys that exhibit IP effects. IP effects are most often approximated as static or modeled in the frequency domain, using frequency-dependent electrical conductivity. It is difficult to translate the frequency-dependent conductivity approach directly to the time domain in a computationally efficient manner. We take an alternative approach in which we model IP relaxations in time using the stretched exponential (SE) function. We incorporate this IP model into a direct time-stepping discretization of the quasistatic time-domain Maxwell equations. We found that modeling of IP effects with this SE approach is asymptotically equivalent to the commonly used Cole-Cole model of IP transformed to the time domain. We have implemented our algorithm using efficient numerical methods that allow it to tackle large-scale problems and are amenable to use in inversion. In particular, we have developed a parallel time-stepping technique that allows us to compute transient electric fields at multiple time steps simultaneously. We demonstrate the behavior of the SE model of IP decay and the efficiency of our algorithm by applying it to synthetic numerical examples that simulate a grounded source IP survey with significant EM effects and a concentric-loop airborne EM sounding over a chargeable body.