Численное решение системы ДУ методом Рунге-Кутты

//Максимальный элемент в массиве
        double max(double[] x)
        {
            double max = x[0];
            for (int i = 1; i < x.GetLength(0); i++)
                if (x[i] > max)
                    max = x[i];
 
            return max;
        }
 
        //l-норма матрицы
        double lnorm(int[][] x)
        {
            int s = 0;
            for (int j = 0; j < 2; j++)
            {
                int max = x[0][j];
                for (int i = 1; i < x.GetLength(0); i++)
                    max = Math.Max(max, x[i][j]);
                s += max;
            }
            return s;
        }
 
        //вектор-функция F, не зависит от времени, поэтому данный параметр не нужен
        int[] f(int[] x, double[] r, double[] k, double[] a)
        {
            int[] result = new int[2];
            double[] precomp = new double[2];
 
            precomp[0] = x[0] * (r[0] - k[0] * x[0] - a[0] * x[1]);
            precomp[1] = x[1] * (r[1] - k[1] * x[0] - a[1] * x[1]);
 
 
            //Проверка получившихся значений: в случае выхода за пределы программа выкинет ошибку, а отрицательным результат быть не может из смысла задачи
            if (precomp[0] > Int32.MaxValue)
                result[0] = Int32.MaxValue;
            else if (precomp[0] < 0)
                result[0] = 0;
            else
                result[0] = Convert.ToInt32(precomp[0]);
 
            if (precomp[1] > Int32.MaxValue)
                result[1] = Int32.MaxValue;
            else if (precomp[1] < 0)
                result[1] = 0;
            else
                result[1] = Convert.ToInt32(precomp[1]);
 
            return result;
        }
 
        //Заметим, что в методе Рунге-Кутта первое приближение k1 совпадает с исходной функцией
 
        //второе приближение: f(t + h/2, y + h/2 * k1);
        int[] k2(int[] x, double[] r, double[] k, double[] a, double h)
        {
            int[] res = new int[2];
            res[0] = Convert.ToInt32(x[0] + h / 2 * f(x, r, k, a)[0]);
            res[1] = Convert.ToInt32(x[1] + h / 2 * f(x, r, k, a)[1]);
 
            return f(res, r, k, a);
 
        }
        
        //третье приближение: f(t + h/2, y + h/2 * k2);
        int[] k3(int[] x, double[] r, double[] k, double[] a, double h)
        {
            int[] res = new int[2];
            res[0] = Convert.ToInt32(x[0] + h / 2 * k2(x, r, k, a, h)[0]);
            res[1] = Convert.ToInt32(x[1] + h / 2 * k2(x, r, k, a, h)[1]);
 
            return f(res, r, k, a);
        }
 
        //четвертое приближение: f(t + h, y + h * k3);
        int[] k4(int[] x, double[] r, double[] k, double[] a, double h)
        {
            int[] res = new int[2];
            res[0] = Convert.ToInt32(x[0] + h * k3(x, r, k, a, h)[0]);
            res[1] = Convert.ToInt32(x[1] + h * k3(x, r, k, a, h)[1]);
 
            return f(res, r, k, a);
        }
 
        //Подсчет на i шаге
        //y_i+1 = y_i + h/6 * (k1 + 2k2 + 2k3 + k4)
        int[] rk4i(int[] x, double[] r, double[] k, double[] a, double h)
        {
            int[] resf = f(x, r, k, a);
            int[] resk2 = k2(x, r, k, a, h);
            int[] resk3 = k3(x, r, k, a, h);
            int[] resk4 = k4(x, r, k, a, h);
            int[] res = new int[2];
 
            res[0] = Convert.ToInt32(x[0] + h / 6 * (resf[0] + 2 * resk2[0] + 2 * resk3[0] + resk4[0]));
            res[1] = Convert.ToInt32(x[1] + h / 6 * (resf[1] + 2 * resk2[1] + 2 * resk3[1] + resk4[1]));
 
            return res;
        }
             
        //Подсчет с шагом h
        int[][] rk4h(int[] x, double[] r, double[] k, double[] a, double h, int T)
        {
            int N = Convert.ToInt32(T / h);
            int[][] res = new int[N][];
            for (int i = 0; i < N; i++)
                res[i] = new int[2];
 
            res[0] = rk4i(x, r, k, a, h);
            for (int i = 1; i < N; i++)
                res[i] = rk4i(res[i - 1], r, k, a, h);
 
            return res;
        }
 
        //Подсчет с точностью eps по правилу Рунге
        int[][] rk4(int[] x, double[] r, double[] k, double[] a, int T, double eps)
        {
            double h1 = 2.5 / max(r);
            double h2 = h1 / 2;
 
            int[][] y1 = rk4h(x, r, k, a, h1, T);
            int[][] y2 = rk4h(x, r, k, a, h2, T);
 
            int N1 = Convert.ToInt32(T / h1);
            int[][] re = new int[N1][];
            for (int i = 0; i < N1; i++)
                re[i] = new int[2];
            for (int i = 0; i < N1; i++)
                for (int j = 0; j < 2; j++)
                    re[i][j] = Math.Abs(y1[i][j] - y2[2 * i][j]);
            double rung = lnorm(re) / 15;                                               //4 порядок точности
 
            while (rung > eps)
            {
                h1 = h2;
                h2 = h1 / 2;
                y1 = rk4h(x, r, k, a, h1, T);
                y2 = rk4h(x, r, k, a, h2, T);
 
                N1 = Convert.ToInt32(T / h1);
                Array.Resize(ref re, N1);
                for (int i = 0; i < N1; i++)
                    re[i] = new int[2];
                for (int i = 0; i < N1; i++)
                    for (int j = 0; j < 2; j++)
                        re[i][j] = Math.Abs(y1[i][j] - y2[2 * i][j]);
                rung = lnorm(re) / 15;
            }
 
            return y1;
        }