A fast numerical method is presented for the simulation of complicated 3D structures, such as inductors constructed from litz or stranded wires, above or sandwiched between planar lossy magnetic media. Basing upon its smoothness, the quasistatic multilayer Green's function is numerically computed using finite differences, and its source height dependence is computed using an optimal Toeplitz-plus-Hankel decomposition. We show that a modified precorrected FFT method can be applied to reduce the dense linear algebra problem to near-linear time, and that frequency-dependent setups can be avoided to result in a considerable speed-up. Experimental verifications are made for a 16-strand litz wire coil, realistically modeled down to each individual strand. Results are obtained in 2-3 hours, showing an excellent agreement to measurements, and can be used to study the impact of transposition patterns in litz wire construction.