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.