Since the beginning of the COVID-19 epidemic there has been an intensive global research effort devoted to the study of the virus and in particular to the development of mathematical models of COVID-19 that involve systems of ordinary differential equations (ODEs). In this paper, we consider systems of ODEs rising from an SEIR epidemic model with interventions. The impact of these interventions is that the solution to the model is nonsmooth at the points in time where the interventions are introduced or removed. This problem is sufficiently challenging that standard ODE solvers are not able to obtain numerical solutions of these models that have even moderate accuracy. However, we show in this paper that dramatic improvements in the accuracy and reliability of approximate solutions of this model can be obtained by employing carefully chosen, robust, numerical ODE methods. In particular, we consider an algorithm that can automatically detect, and efficiently and accurately handle the discontinuities that arise when the model includes interventions that are imposed in an attempt to restrict the spread of the virus and later removed when the spread of the virus is diminished.