  We present an asymptotically and unconditionally stable numerical method to
account for the momentum transfer between multiple species. Momentum is
conserved to machine precision. This implies that the asymptotic equilibrium
corresponds to the velocity of the center of mass. Aimed at studying dust
dynamics, we implement this numerical method in the publicly available code
FARGO3D. To validate our implementation, we develop a test suite for an
arbitrary number of species, based on analytical or exact solutions of problems
related to perfect damping, damped sound waves, shocks, local and global
gas-dust radial drift in a disk and, linear streaming instability. In
particular, we obtain first-order, steady-state solutions for the radial drift
of multiple dust species in protoplanetary disks, in which the pressure
gradient is not necessarily small. We additionally present non-linear
shearing-box simulations of the streaming instability and compare them with
previous results obtained with Lagrangian particles. We successfully validate
our implementation by recovering the solutions from the test suite to second-
and first-order accuracy in space and time, respectively. From this, we
conclude that our scheme is suitable, and very robust, to study the
self-consistent dynamics of several fluids. In particular, it can be used for
solving the collisions between gas and dust in protoplanetary disks, with any
degree of coupling.
