differentiating the SEIR model in Java

I'm attempting to make a simulation of the SEIR epidemic model.
It contains four parts:

• Susceptibles (non-infected)
• Exposed (infected but not infectious yet)
• Infectious (infected and infectious)
• Removed (recovered/dead)

where gamma γ is the infection rate and beta β is the reovery/death rate.

I've previously used the SIR model, a more basic model where E and I are combined, which uses these equations:

From another thread I've used a solution to simulate SIR using this code:

double dS = (beta * S.get(day) * I.get(day) / N);
double newS = (S.get(day) - dS);
double newI = (I.get(day) + dS - gamma * I.get(day));
double newR = (R.get(day) + gamma * I.get(day));

This works fine using the Euler's method. However, I've tried to manipulate this to try and fit the SEIR model (which has these equations:)

where u is the death rate, delta is the birth rate and a is the incubation period. I've made an attempt to try and use a similar method to work for SEIR but I'm unsuccessful to simulate it well at all. This isn't really a problem with the variables but as a whole differentiating these complex equations. Wondering if anyone could help, thanks.

Really should've realised this earlier but from messing around with random sign changes, worked out that everything apart from 'newS' requires getting the previous day's number and plussing the new dS, rather than minusing it. My SIR code already did this. Don't really know how I missed this..

New working code:

int totalDays = 160; // How many days/times to loop
int N = 1000; // Population
int I0 = 1; // Starting infected/exposed
double beta = 0.2; // Infection rate
double gamma = 1.0/10.0; // recovery time (days to the -1)
double a = 1.0/2.0; // incubation period (days to the -1)
List<Double> S = new ArrayList<>();
List<Double> E = new ArrayList<>();
List<Double> I = new ArrayList<>();
List<Double> R = new ArrayList<>();
private void createData() {
final int R0 = 0;
final int S0 = N - E0 - R0;
S.add((double) S0);
E.add((double) I0);
I.add(0.0);
R.add(0.0);
for (int day = 1; day < totalDays + 1; day++) {
double[] derivative = deriv(day);
S.add(derivative[0]);
E.add(derivative[1]);
I.add(derivative[2]);
R.add(derivative[3]);
}
}
private double[] deriv(int day) {
day = day - 1;
double dS = (beta * S.get(day) * I.get(day)) / N;
double newS = S.get(day) - (dS);
double newE = E.get(day) + (dS - (a * E.get(day)));
double newI = I.get(day) + ((a * E.get(day)) - (gamma * I.get(day)));
double newR = R.get(day) + (gamma * I.get(day));
return new double[] {newS, newE, newI, newR};
}