Рунге-Кутта (RK4) для системы дифференциальных уравнений на Java

Это предложение в основном является результатом этой темы: Дифференциальные уравнения в Java.
В основном, Я попытался последовать совету Джейсона С. и реализовать численные решения дифференциальных уравнений с помощью метода Рунге-Кутта (RK4).

Всем привет, пытаюсь создать простую имитационную программу модели SIR-эпидемий на java. По сути, SIR определяется системой трех дифференциальных уравнений:
S '(t) = - lamda (t) * S (t)
I' (t) = lamda (t) * S (t) - гамма (t) * I (t)
R '(t) = гамма (t) * I (t)
S - восприимчивые люди, I - инфицированные люди, R - выздоровевшие люди. lamda (t) = [c * x * I (t)] / N (T) c - количество контактов, x - инфекционность (вероятность заболеть после контакта с больным человеком), N (t) - общая популяция (которая постоянна).
гамма (t) = 1 / продолжительность болезни (константа)

После первой не очень успешной попытки я попытался решить эти уравнения с помощью Runge-KUtta, и эта попытка привела к следующему коду:

package test;

public class Main {
    public static void main(String[] args) {


        double[] S = new double[N+1];
        double[] I = new double[N+1];
        double[] R = new double[N+1];

        S[0] = 99;
        I[0] = 1;
        R[0] = 0;

        int steps = 100;
        double h = 1.0 / steps;
        double k1, k2, k3, k4;
        double x, y;
        double m, n;
        double k, l;

        for (int i = 0; i < 100; i++) {
            y = 0;
            for (int j = 0; j < steps; j++) {
                x = j * h;
                k1 = h * dSdt(x, y, S[j], I[j]);
                k2 = h * dSdt(x+h/2, y +k1/2, S[j], I[j]);
                k3 = h * dSdt(x+h/2, y+k2/2, S[j], I[j]);
                k4 = h * dSdt(x+h, y + k3, S[j], I[j]);
                y += k1/6+k2/3+k3/3+k4/6;
            }
            S[i+1] = S[i] + y;
            n = 0;
            for (int j = 0; j < steps; j++) {
                m = j * h;
                k1 = h * dIdt(m, n, S[j], I[j]);
                k2 = h * dIdt(m+h/2, n +k1/2, S[j], I[j]);
                k3 = h * dIdt(m+h/2, n+k2/2, S[j], I[j]);
                k4 = h * dIdt(m+h, n + k3, S[j], I[j]);
                n += k1/6+k2/3+k3/3+k4/6;
            }
            I[i+1] = I[0] + n;
            l = 0;
            for (int j = 0; j < steps; j++) {
                k = j * h;
                k1 = h * dRdt(k, l, I[j]);
                k2 = h * dRdt(k+h/2, l +k1/2, I[j]);
                k3 = h * dRdt(k+h/2, l+k2/2, I[j]);
                k4 = h * dRdt(k+h, l + k3, I[j]);
                l += k1/6+k2/3+k3/3+k4/6;
            }
            R[i+1] = R[i] + l;
        }
        for (int i = 0; i < 100; i ++) {
            System.out.println(S[i] + " " + I[i] + " " + R[i]);
        }
    }

    public static double dSdt(double x, double y, double s, double i) {
        return (- c * x * i / N) * s;
    }
    public static double dIdt(double x, double y, double s, double i) {
        return (c * x * i / N) * s - g * i;
    }
    public static double dRdt(double x, double y, double i) {
        return g*i;
    }

    private static int N = 100;

    private static int c = 5;
    private static double x = 0.5;      
    private static double g = (double) 1 / x;
}

Это не работает, потому что количество больных (I) должно сначала увеличиваться, а затем уменьшаться примерно до 0, а число выздоровевших должно строго увеличиваться. Общее количество больных + здоровых + выздоровевших должно быть 100, но мой код дает некоторые странные результаты:

99.0 1.0 0.0  
98.9997525 0.9802475 0.03960495  
98.99877716805084 0.9613703819491656 0.09843730763898331  
98.99661215494893 0.9433326554629141 0.1761363183872249  
98.99281287394908 0.9261002702516101 0.2723573345404987  
98.98695085435723 0.9096410034385773 0.3867711707625441  
98.97861266355956 0.8939243545756241 0.5190634940761019  
98.96739889250752 0.8789214477384787 0.6689342463444292  
98.95292320009872 0.8646049401404658 0.8360970974155659  
98.93481141227473 0.8509489367528628 1.0202789272217598  
98.91270067200323 0.8379289104653137 1.22121933523726  
98.8862386366277 0.8255216273600343 1.438670175799961
98.8550827193552 0.8137050767097959 1.672395117896858  

Не могу найти ошибку, посоветуйте пожалуйста! Спасибо заранее!