Решить приближенно задачу Коши для ОДУ 3 порядка

using System;
using System.Collections.Generic;
 
namespace rwsh_rk4
{
    abstract class TRungeKutta 
    {
        public int N;
        double t; // текущее время 
        public double[] Y; // искомое решение Y[0] - само решение, Y[i] - i-тая производная решения
 
        double[] YY, Y1, Y2, Y3, Y4; // внутренние переменные 
 
        public TRungeKutta(int N) // N - размерность системы 
        {
            this.N = N; // сохранить размерность системы
 
            if (N < 1)
            {
                this.N = -1; // если размерность меньше единицы, то установить флаг ошибки
                return; // и выйти из конструктора
            }
 
            Y = new double[N]; // создать вектор решения
            YY = new double[N]; // и внутренних решений
            Y1 = new double[N];
            Y2 = new double[N];
            Y3 = new double[N];
            Y4 = new double[N];
        }
 
        public void SetInit(double t0, double[] Y0) // установить начальные условия.
        {                                           // t0 - начальное время, Y0 - начальное условие
            t = t0;
            int i;
            for (i = 0; i < N; i++)
            {
                Y[i] = Y0[i];
            }
        }
 
        public double GetCurrent() // вернуть текущее время
        {
            return t;
        }
 
        abstract public void F(double t, double[] Y, ref double[] FY); // правые части системы.
 
        public void NextStep(double dt) // следующий шаг метода Рунге-Кутта, dt - шаг по времени (может быть переменным)
        {
            if(dt<0)
            {
                return;
            }
 
            int i;
 
            F(t, Y, ref Y1); // рассчитать Y1
 
            for (i = 0; i < N; i++)
            {
                YY[i] = Y[i] + Y1[i] * (dt / 2.0);
            }
            F(t + dt / 2.0, YY, ref Y2); // рассчитать Y2
 
            for (i = 0; i < N; i++)
            {
                YY[i] = Y[i] + Y2[i] * (dt / 2.0);
            }
            F(t + dt / 2.0, YY, ref Y3); // рассчитать Y3
 
            for (i = 0; i < N; i++)
            {
                YY[i] = Y[i] + Y3[i] * dt;
            }
            F(t + dt, YY, ref Y4); // рассчитать Y4
 
            for (i = 0; i < N; i++)
            {
                Y[i] = Y[i] + dt / 6.0 * (Y1[i] + 2.0 * Y2[i] + 2.0 * Y3[i] + Y4[i]); // рассчитать решение на новом шаге
            }
 
            t = t + dt; // увеличить шаг
 
        }
    }
 
    class TMyRK : TRungeKutta
    {
        public TMyRK(int aN) : base(aN) { }
 
        public override void F(double t, double[] Y, ref double[] FY) 
        {
            FY[0] = Y[1]; // пример математический маятник 
            FY[1] = -Y[0]; // y''(t)+y(t)=0
        }
    }
 
    class Program
    {
        static void Main(string[] args)
        {
            TMyRK RK4 = new TMyRK(2);
 
            double[] Y0 = {0, 1}; // зададим начальные условия y(0)=0, y'(0)=1
 
            RK4.SetInit(0, Y0);
 
            while (RK4.GetCurrent() < 10) // решаем до 10
            {
                Console.WriteLine("{0}\t{1}\t{2}", RK4.GetCurrent(), RK4.Y[0], RK4.Y[1]); // вывести t, y, y'
 
                RK4.NextStep(0.01); // расчитать на следующем шаге, шаг интегрирования dt=0.01
            }
        }
    }
}