Runge-kutta 4th, SIR epidemic model - C#

Формулировка задачи:

Есть код:

namespace MMHTP
{
  public partial class Form1 : Form
  {
    public Form1()
    {
      InitializeComponent();
    }
    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;
 
    private void button1_Click(object sender, EventArgs e)
    {
      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 axisX = h;
      double k1, k2, k3, k4;
      double x, y;
      double m, n;
      double k, l;
      chart1.Series[0].ChartType = System.Windows.Forms.DataVisualization.Charting.SeriesChartType.Spline;
      chart1.Series[0].Points.Clear();
 
      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;
        //chart1.Series[0].Points.AddXY(axisX, S[i+1]);
        //axisX += h;
      }
      for (int i = 0; i < 100; i++)
      {
        Console.WriteLine(S[i] + " " + I[i] + " " + R[i]);
        Console.WriteLine(S[i] + I[i] + R[i]);
        // тут сами поизменяет/повыбирайте тип вывода графика
 
      }
      Console.WriteLine("end");
    }
  }
}

Должно соблюдаться тождество S'+I'+R' = 0 но на выходе получаются какие-то странные значения:

99 1 0 98,9997525 0,9802475 0,03960495 98,9987771680508 0,961370381949166 0,0984373076389833 98,9966121549489 0,943332655462914 0,176136318387225 98,9928128739491 0,92610027025161 0,272357334540499 98,9869508543572 0,909641003438577 0,386771170762544 98,9786126635596 0,893924354575624 0,519063494076102 98,9673988925075 0,878921447738479 0,668934246344429 98,9529232000987 0,864604940140466 0,836097097415566 98,9348114122747 0,850948936752863 1,02027892722176 98,9127006720032 0,837928910465314 1,22121933523726 98,8862386366277 0,825521627360034 1,43867017579996 98,8550827193552 0,813705076709796 1,67239511789686 98,818899371916 0,802458405342307 1,9221692281006 98,7773634056693 0,791761856042978 2,1877785754252 98,7301573486501 0,781596709694539 2,4690198569437 98,6769708362558 0,771945230875826 2,76570004307971 98,6175000334562 0,762790616663641 3,07763604154899 98,5514470865859 0,754116948400998 3,40465437898629 98,4785196029359 0,745909146212691 3,74659089934785 98,3984301565093 0,738152926064976 4,1032904782307 98,3108958184459 0,73083475918055 4,47460675229717 98,2156377107465 0,723941833632996 4,8604018630363 98,1123805820507 0,717462017956608 5,26054621413456 98,0008524043343 0,711383826618149 5,67491824176518 97,880783989497 0,705696387206703 6,10340419713994 97,7519086249158 0,700389409206458 6,54589794069882 97,6139617271318 0,695453154225135 7,00230074734221 97,4666805129306 0,690878407557819 7,47252112213675 97,3098036871636 0,6866564509724 7,95647462595075 97,1430711467411 0,682779036608538 8,45408371049691 96,9662237003076 0,679238361887309 8,96527756228082 96,7790028031904 0,676027045333337 9,48999195497139 96,5811503072844 0,673138103215414 10,0281691097263 96,3724082256141 0,670564926915411 10,5797575630195 96,1525185113823 0,668301260938639 11,1447120415314 95,9212228513884 0,666341181481902 11,722993343673 95,6782624737686 0,664679075478191 12,3145682273242 95,423377970078 0,663309620039456 12,9194093033762 95,1563091318049 0,662227762221128 13,5374949346725 94,8767948014745 0,661428699034077 14,1688091399496 94,5845727385658 0,660907857631593 14,8133415023793 94,2793795005354 0,660660875600693 15,471087082321 93,960950339305 0,660683581288711 16,1420463338885 93,6290191136398 0,66097197409773 16,8262250249379 93,2833182179094 0,661522204680969 17,523634160081 92,9235785277895 0,662330554976848 18,2342899063236 92,5495293635288 0,663393418018105 18,9582135209265 92,1608984714707 0,664707277455097 19,6954312810786 91,7574120245843 0,666268686734361 20,4459744149653 91,338794642822 0,668074247875595 21,2098790338096 90,9047694341851 0,670120589792632 21,9871860644497 90,4550580574384 0,672404346106641 22,7779411820119 89,9893808074733 0,67492213240286 23,5821947422222 89,5074567243773 0,677670522885644 24,4000017128902 89,0090037273187 0,680646026390587 25,231421604086 88,4937387744059 0,683845061717021 26,0765183965162 87,9613780497263 0,687263932249366 26,9353604675914 87,4116371788095 0,690898799841681 27,8080205146633 86,8442314737911 0,694745657946415 28,6945754748942 86,2588762095843 0,698800303975884 29,5951064412047 85,6552869323805 0,703058310893419 30,509698573733 85,0331798018116 0,707514998040605 31,4384410062221 84,3922719681049 0,712165401217563 32,3814267467355 83,7322819855465 0,717004242044912 33,3387525720899 83,0529302635432 0,722025896648979 34,3105189153772 82,3539395565298 0,727224363726036 35,2968297459391 81,6350354939111 0,73259323205688 36,2977924411421 80,8959471511481 0,73812564756002 37,3135176492963 80,1364076630038 0,743814279990073 38,3441191430503 79,356154879842 0,749651289407751 39,3897136625925 78,5549320677309 0,755628292569007 40,450420747986 77,7324886529338 0,761736329403486 41,5263625599676 76,8885810111759 0,767965829776349 42,6176636885447 76,0229733018461 0,774306580752663 43,7244509487369 75,1354383470441 0,780747694609821 44,8468531628213 74,225758555089 0,787277577870664 45,9850009284631 73,2937268877893 0,793883901657906 47,139026372138 72,3391478704155 0,80055357369888 48,309062887287 71,3618386429284 0,807272712338204 49,4952448566826 70,3616300505883 0,814026622944331 50,6977073585372 69,3383677716101 0,820799777123657 51,9165858559343 68,2919134790306 0,827575795182504 53,152015869235 67,2221460334277 0,834337432302156 54,4041326311817 66,1289627025662 0,84106656891478 55,6730707245068 65,0122804034534 0,847744205787679 56,9589637019476 63,8720369616734 0,854350464339223 58,2619436886751 62,7081923822246 0,860864592721224 59,5821409672571 61,5207301254341 0,86726497820854 60,9196835454033 60,3096583808513 0,873529166436556 62,2746967068782 59,0750113313566 0,879633888019866 63,6473025461135 57,816850399051 0,885555093070222 65,0376194872103 56,5352654638467 0,89126799410761 66,4457617881891 55,2303760450435 0,896747117824374 67,8718390315245 53,9023324355902 0,901966366117881 69,3159556021822 52,551316778181 0,906899086751542 70,7782101545749 51,1775440718517 0,911518153936514 72,2586950700464 49,7812630973335 0,915796059046804 73,7574959066988 48,3627572490927 0,919705011588337 75,2746908435829 46,9223452617705 0,923217050438058 76,8103501214759

Не могу разобраться в чем ошибка.

Решение задачи: «Runge-kutta 4th, SIR epidemic model»

textual
                            Листинг программы 
                        

                            public partial class Form1 : Form
{
    public Form1()
    {
        InitializeComponent();
    }
 
    IEnumerable<Tuple<double, Vector>> rgk(double x, Vector y, Func<double, Vector, Vector> F, double h)
    {
        while (true)
        {
            yield return Tuple.Create(x, y);
            var k1 = h * F(x + 0.0 * h, y);
            var k2 = h * F(x + 0.5 * h, y + k1 / 2);
            var k3 = h * F(x + 0.5 * h, y + k2 / 2);
            var k4 = h * F(x + 1.0 * h, y + k3 / 1);
            y += (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            x += h;
        }
    }
 
    private void button1_Click(object sender, EventArgs e)
    {
        chart1.Series[0].Points.Clear();
        chart1.Series[1].Points.Clear();
        chart1.Series[2].Points.Clear();
        chart1.Series[3].Points.Clear();
 
        int N = 100; // Число людей
        double c = 50; // Коэффициент взаимодействия между людьми
        double gx = 0.5;
        double g = 1.0 / gx;
 
        Func<double, Vector, Vector> df = (x, f) => new Vector(
            /* dS = */ ((-c * x * f.I) / N) * f.S,
            /* dI = */ ((c * x * f.I) / N) * f.S - g * f.I,
            /* dR = */ g * f.I
            );
 
        int n = 100; // размер сетки для интегрирования и графика одновременно
        double a = 0, b = 1; // её границы
        double h = (b - a) / n; // шаг
        foreach (var y in rgk(a, new Vector(99, 1, 0), df, h).Take(n + 1))
        {
            chart1.Series[0].Points.AddXY(y.Item1, y.Item2.S);
            chart1.Series[1].Points.AddXY(y.Item1, y.Item2.I);
            chart1.Series[2].Points.AddXY(y.Item1, y.Item2.R);
            chart1.Series[3].Points.AddXY(y.Item1, y.Item2.N);
            Debug.Print("Время: {0,2:F} {1}", y.Item1, y.Item2);
        }
    }
 
    struct Vector
    {
        double[] c;
        public Vector(params double[] Coords) { c = Coords; }
        public static Vector operator +(Vector a, Vector b) => new Vector(a.c.Select((x, i) => x + b.c[i]).ToArray());
        public static Vector operator *(double a, Vector b) => new Vector(b.c.Select((x, i) => a * x).ToArray());
        public static Vector operator /(Vector a, double b) => new Vector(a.c.Select((x, i) => x / b).ToArray());
 
        public double S => c[0];
        public double I => c[1];
        public double R => c[2];
        public double N => S + I + R;
 
        public override string ToString()
        {
            return string.Format(
                "Здоровых: {0,5:F} Больных: {1,5:F} Переболевших: {2,5:F} Всего: {3,5:F}",
                S, I, R, N);
        }
    }
}
                        