  Ambipolar diffusion is a physical mechanism related to the drift between
charged and neutral particles in a partially ionized plasma that is key in many
different astrophysical systems. However, understanding its effects is
challenging due to basic uncertainties concerning relevant microphysical
aspects and the strong constraints it imposes on the numerical modeling. Our
aim is to introduce a numerical tool that allows us to address complex problems
involving ambipolar diffusion in which, additionally, departures from
ionization equilibrium are important or high resolution is needed. The primary
application of this tool is for solar atmosphere calculations, but the methods
and results presented here may also have a potential impact on other
astrophysical systems. We have developed a new module for the stellar
atmosphere Bifrost code that improves its computational capabilities of the
ambipolar diffusion term in the Generalized Ohm's Law. This module includes,
among other things, collision terms adequate to processes in the coolest
regions in the solar chromosphere. As a key feature of the module, we have
implemented the Super Time-Stepping (STS) technique, that allows an important
acceleration of the calculations. We have also introduced hyperdiffusion terms
to guarantee the stability of the code. We show that to have an accurate value
for the ambipolar diffusion coefficient in the solar atmosphere it is necessary
to include as atomic elements in the equation of state not only hydrogen and
helium but also the main electron donors like sodium, silicon and potassium. In
addition, we establish a range of criteria to set up an automatic selection of
the free parameters of the STS method that guarantees the best performance,
optimizing the stability and speed for the ambipolar diffusion calculations. We
validate the STS implementation by comparison with a self-similar analytical
solution.
