We prove the global non-linear stability, without symmetry assumptions, of slowly rotating charged black holes in de Sitter spacetimes in the context of the initial value problem for the Einstein-Maxwell equations: If one perturbs the initial data of a slowly rotating Kerr-Newman-de Sitter (KNdS) black hole, then in a neighborhood of the exterior region of the black hole, the metric and the electromagnetic field decay exponentially fast to their values for a possibly different member of the KNdS family. This is a continuation of recent work of the author with Vasy on the stability of the Kerr-de Sitter family for the Einstein vacuum equations. Our non-linear iteration scheme automatically finds the final black hole parameters as well as the gauge in which the global solution exists; we work in a generalized wave coordinate/Lorenz gauge, with gauge source functions lying in a suitable finite-dimensional space. We include a self-contained proof of the linear mode stability of Reissner-Nordstr\"om-de Sitter black holes, building on work by Kodama-Ishibashi. In the course of our non-linear stability argument, we also obtain the first proof of the linear (mode) stability of slowly rotating KNdS black holes using robust perturbative techniques.