  Many important applications are modelled by differential equations with
positive solutions. However, it remains an outstanding open problem to develop
numerical methods that are both (i) of a high order of accuracy and (ii)
capable of preserving positivity. It is known that the two main families of
numerical methods, Runge-Kutta methods and multistep methods, face an order
barrier: if they preserve positivity, then they are constrained to low
accuracy: they cannot be better than first order. We propose novel methods that
overcome this barrier: our methods are of second order, and they are guaranteed
to preserve positivity. Our methods apply to a large class of differential
equations that have a special graph Laplacian structure, which we elucidate.
The equations need be neither linear nor autonomous and the graph Laplacian
need not be symmetric. This algebraic structure arises naturally in many
important applications where positivity is required. We showcase our new
methods on applications where standard high order methods fail to preserve
positivity, including infectious diseases, Markov processes, master equations
and chemical reactions.
