Charged dust in magnetized gas is subject to aerodynamic drag and the Lorentz force, both of which can be difficult to model numerically in a fluid simulation where their characteristic timescales are shorter than the timestep set by the Courant condition. We present a two-stage general implicit Runge-Kutta (GIRK) source integrator for multifluid charged dust coupled to an ideal magnetohydrodynamic (MHD) gas, implemented in the open source GPU-accelerated quokka code. Drag and Lorentz forces are advanced together in one implicit source update, avoiding errors in the terminal drift that can arise when the two forces are advanced in separate substeps. The timestep is compared with a drag-Lorentz timescale to select between resolved and stiff branches. The resolved branch uses Gauss-Legendre coefficients, while the stiff branch uses a coefficient set designed to recover the terminal drift set by the drag-Lorentz balance. We also allow the dust stopping time and charge-to-mass ratio to depend on the evolving state, evaluating them separately at the two implicit stage states and solving the resulting nonlinear stage equations by Picard iteration. The source update conserves total gas-dust momentum, and its computational cost scales linearly with the number of dust species. Tests reproduce the expected drag damping and terminal drift, confirm the predicted gyromotion amplitude and phase errors, demonstrate accurate treatment of state-dependent source coefficients, and recover the expected wave, shock, and instability dynamics in dusty MHD flows.