  We describe a method for incorporating ambipolar diffusion in the strong
coupling approximation into a multidimensional magnetohydrodynamics code based
on the total variation diminishing scheme. Contributions from ambipolar
diffusion terms are included by explicit finite difference operators in a fully
unsplit way, maintaining second order accuracy. The divergence-free condition
of magnetic fields is exactly ensured at all times by a flux-interpolated
constrained transport scheme. The super time stepping method is used to
accelerate the timestep in high resolution calculations and/or in strong
ambipolar diffusion. We perform two test problems, the steady-state oblique
C-type shocks and the decay of Alfv\'en waves, confirming the accuracy and
robustness of our numerical approach. Results from the simulations of the
compressible MHD turbulence with ambipolar diffusion show the flexibility of
our method as well as its ability to follow complex MHD flows in the presence
of ambipolar diffusion. These simulations show that the dissipation rate of MHD
turbulence is strongly affected by the strength of ambipolar diffusion.
