Time-domain problems for controlled-source electromagnetic exploration require accurate discretization of the solution for multiple spacial and temporal scales. Therefore, forward simulation using conventional computational methods becomes computationally expensive, even without accounting for induced-polarization (IP) effects. These effects create another complication caused by the presence of a convolution integral in the time-domain Maxwell system. We suggested a novel, fast, and robust algorithm to solve the 3D time-domain electromagnetic (EM) problems that can be considered as a generalization of the spectral Lanczos decomposition method. The new method also allowed us to incorporate the IP effects without significant cost increase. The discretized large-scale Maxwell system was projected onto a small subspace consisting of the Laplace-domain solutions (the so-called parameter-dependent Krylov subspace) for an optimally chosen set of Laplace parameters. The projected system preserved stability and passivity of the original problem. Moreover, our approach (even without the IP effects) yielded an optimal solution within a wide class of computational algorithms that included the conventional time-domain finite-difference, discrete Fourier transform and spectral Lanczos decomposition methods. Numerical examples for the controlled-source EM problem showed that the new algorithm produces accurate solutions on time intervals spanning from milliseconds to hundreds of seconds with the cost of (at most) 60 time steps of the implicit finite-difference time domain scheme. This showed significant improvement even compared with results for nonpolarized media reported in recent literature. Additionally, the new algorithm had the unique capability to accurately handle large-scale 3D models, including the IP effects.