  Numerical simulations of dust-gas dynamics are one of the fundamental tools
in astrophysical research, such as the study of star and planet formation. It
is common to find tightly coupled dust and gas in astrophysical systems, which
demands that any practical integration method be able to take time steps
$\Delta t$ much longer than the stopping time $t_{\rm s}$ due to drag. A number
of methods have been developed to ensure stability in this stiff ($\Delta t\gg
t_{\rm s}$) regime, but there remains large room for improvement in terms of
accuracy. In this paper, we describe an easy-to-implement method, the
"staggered semi-analytic method" (SSA), and conduct numerical tests to compare
it to other implicit and semi-analytic methods, including the $2^{\rm nd}$
order implicit method and the Verlet method. SSA makes use of a staggered step
to better approximate the terminal velocity in the stiff regime. In
applications to protoplanetary disks, this not only leads to
orders-of-magnitude higher accuracy than the other methods, but also provides
greater stability, making it possible to take time steps 100 times larger in
some situations. SSA is also $2^{\rm nd}$ order accurate and symplectic when
$\Delta t \ll t_{\rm s}$. More generally, the robustness of SSA makes it
applicable to linear dust-gas drag in virtually any context.
