SEIR equations solver

This webpage uses the Runge-Kutta-Fehlberg fourth-order method with fifth-order error checking (RKF45) to approximate the solution to the SEIR equations with the $\delta$ parameter to account for quarantine effects: $$\begin{aligned} \frac{dS}{dt} & = \Lambda N - \mu S - \frac{\beta I (1-\delta)S}{N} \\ \frac{dE}{dt} & = \frac{\beta I (1-\delta) S}{N} - (\mu +a ) E \\ \frac{dI}{dt} & = a E - (\gamma +\mu ) I \\ \frac{dR}{dt} & = \gamma I - \mu R. \end{aligned}$$ Where $S$ is the number of susceptible persons, $E$ is the number of exposed persons, $I$ is the number of infectious persons and $R$ is the number of recovered persons. $a$ is the inverse of the average incubation period. $\beta$ is a parameter that pertains to the average number of contacts per person per time and the rate of transmission for the disease. $\gamma$ is the inverse of the average time a person is infected with the disease. $\Lambda$ is the birth rate. $\mu$ is the overall population death rate (not only including the disease death rate). $N$ is the total population. My original model had $\gamma I$ multiplied by $1-\delta$, but as quarantine should not affect how long it takes for people to recover, it should not affect this term.