Late-time negative responses in central-loop transient electromagnetic (TEM) data are often linked to the induced polarization (IP) effect. Early methods for modeling the IP effect in TEM data try to avoid calculating the fractional derivative arising from considering the Cole-Cole model by either using the Fourier transform to convert frequency-domain responses to the time domain or approximating the fractional derivative in the time domain directly. The frequency-to-time conversion method suffer from accuracy issues if the number of frequencies calculated is small. The time-domain approximation method also has accuracy issues because of simplified Cole-Cole models. The Caputo series can approximate fractional derivatives accurately if historic electromagnetic (EM) fields are saved.However, the storage of historic EM fields leads to a significant memory consumption. We introduce the sum-of-exponentials (SOE) method to discretize fractional derivatives, which does not need to store field values from previous times except for the first two time-steps. We discretize the resulting partial differential equations from the SOE discretization using a finite-difference time-domain (FDTD) approach. Additionally, we improve computational efficiency by employing the direct-splitting strategy to transform large sparse matrices into smaller diagonally dominant tridiagonal matrices. We validate the accuracy and efficiency of our algorithm by comparing it with the Caputo approximation method using a chargeable half-space model. Furthermore, we compare our results with existing literature data for a chargeable anomaly in a nonchargeable half space. Finally, we analyze the response characteristics of the IP effect using a block-in-half space model.