Geophysical exploration methods that use controlled electromagnetic sources in time-domain are becoming ubiquitous because of their ease of deployment and data coverage capabilities. While these prospection methodologies have evolved significantly, most computational numerical modeling techniques have been developed in the frequency domain. Given this evolution, it is pertinent to ask whether some known advantages of implicit finite-difference modeling techniques, that are common among other disciplines, apply to the TDEM problem in geophysics. In order to explore the potential advantages of 3D time-domain modeling for TDEM we analyze the differences between two implicit finite-difference formulations: a single order Backward Euler scheme and a 2nd order Crank-Nicolson scheme. To validate our algorithms we tested them with existing analytical solutions for simple geometries for various electromagnetic sources. The results show an acceptable match of both numerical schemes to the analytical solutions and higher accuracy of the 2nd order discrete operator without significant difference in computing time.