We consider a time-dependent diffusion-reaction model for two vector unknowns, satisfying a divergence-free constraint, and the associated scalar Lagrange multiplier. The motivation for studying such a model is provided by a plasma physics problem arising in the modeling of nuclear fusion devices (Braginskii equations), where the two vector unknowns represent ion and electron velocities, the scalar unknown is the electrostatic potential and the divergence-free constraint reflects the physical assumption of quasi-neutrality. We first recast the problem in a form reminiscent of the standard Stokes problem, which allows us to recognize the importance of using a compatible discretization for the vector and scalar unknowns, then propose and analyze a stable finite element formulation. Following this, we address some peculiar geometrical aspects of the model, showing how they can be naturally dealt with within our formulation, and finally discuss a solution procedure for the resulting linear system based on the classical Uzawa algorithm. Some numerical experiments complete the paper.