  Drift-diffusion models that account for the motion of both electronic and
ionic charges are important tools for explaining the hysteretic behaviour and
guiding the development of metal halide perovskite solar cells. Furnishing
numerical solutions to such models for realistic operating conditions is
challenging owing to the extreme values of some of the parameters. In
particular, those characterising (i) the short Debye lengths (giving rise to
rapid changes in the solutions across narrow layers), (ii) the relatively large
potential differences across devices and (iii) the disparity in timescales
between the motion of the electronic and ionic species give rise to significant
stiffness. We present a finite difference scheme with an adaptive time step
that is posed on a non-uniform staggered grid that provides second order
accuracy in the mesh spacing. The method is able to cope with the stiffness of
the system for realistic parameters values whilst providing high accuracy and
maintaining modest computational costs. For example, a transient sweep of a
current-voltage curve can be computed in only a few minutes on a standard
desktop computer.
