Уравнение методом Рунге-Кутты - C#

namespace RungeKutta
{
    class Program
    {
        public delegate double Function(double x, double y);
 
        #region Ordinary Differential Equations - Methods
 
        public static double ODE_RungeKutta2(Function f, double x0, double y0, double h, double x)
        {
            double xnew, ynew, k1, k2, result = double.NaN;
            if (x == x0)
                result = y0;
            else if (x > x0)
            {
                do
                {
                    if (h > x - x0) h = x - x0;
                    k1 = h * f(x0, y0);
                    k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
                    ynew = y0 + k2;
                    xnew = x0 + h;
                    x0 = xnew;
                    y0 = ynew;
                } while (x0 < x);
                result = ynew;
            }
            return result;
        }
 
        public static double ODE_RungeKutta4(Function f, double x0, double y0, double h, double x)
        {
            double xnew, ynew, k1, k2, k3, k4, result = double.NaN;
            if (x == x0)
                result = y0;
            else if (x > x0)
            {
                do
                {
                    if (h > x - x0) h = x - x0;
                    k1 = h * f(x0, y0);
                    k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
                    k3 = h * f(x0 + 0.5 * h, y0 + 0.5 * k2);
                    k4 = h * f(x0 + h, y0 + k3);
                    ynew = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
                    xnew = x0 + h;
                    x0 = xnew;
                    y0 = ynew;
                } while (x0 < x);
                result = ynew;
            }
            return result;
        }
 
        public static double ODE_RungeKuttaFehlberg(Function f, double x0, double y0, double x, double h, double tolerance)
        {
            double xnew, ynew, hnew, k1, k2, k3, k4, k5, k6;
            double hmin = 0.0001;
            double hmax = 0.5;
            if (h > hmax) h = hmax;
            if (h < hmin) h = hmin;
 
            while (x0 < x)
            {
                k1 = h * f(x0, y0);
                k2 = h * f(x0 + 0.25 * h, y0 + 0.25 * k1);
                k3 = h * f(x0 + 3 * h / 8, y0 + 3 * k1 / 32 + 9 * k2 / 32);
                k4 = h * f(x0 + 12 * h / 13, y0 + 1932 * k1 / 2197 - 7200 * k2 / 2197 + 7296 * k3 / 2197);
                k5 = h * f(x0 + h, y0 + 439 * k1 / 216 - 8 * k2 + 3680 * k3 / 513 - 845 * k4 / 4104);
                k6 = h * f(x0 + 0.5 * h, y0 - 8 * k1 / 27 + 2 * k2 - 3544 * k3 / 2565 + 1859 * k4 / 4104 - 11 * k5 / 40);
                double error = Math.Abs(k1 / 360 - 128 * k3 / 4275 - 2197 * k4 / 75240 + k5 / 50 + 2 * k6 / 55) / h;
                double s = Math.Pow(0.5 * tolerance / error, 0.25);
                if (error < tolerance)
                {
                    ynew = y0 + 25 * k1 / 216 + 1408 * k3 / 2565 + 2197 * k4 / 4104 - 0.2 * k5;
                    xnew = x0 + h;
                    x0 = xnew;
                    y0 = ynew;
                }
                if (s < 0.1) s = 0.1;
                if (s > 4) s = 4;
                hnew = h * s;
                h = hnew;
                if (h > hmax) h = hmax;
                if (h < hmin) h = hmin;
                if (h > x - x0) h = x - x0;
            } return y0;
        }
        #endregion
 
        #region Ordinary Differential Equations - Test functions
 
        static double f(double x, double y)
        {
            return y * Math.Cos(x);
        }
        #endregion   
 
        #region Ordinary Differential Equations - Test Methodes
 
        static void TestRungeKutta2()
        {
            double h = 0.001;
            double x0 = 0.0;
            double y0 = 1.0;
            Console.WriteLine("\n Results from the 2nd-order Runge-Kutta method with h = {0}\n", h);
            double result = y0;
            for (int i = 0; i < 11; i++)
            {
                double x = 0.1 * i;
                result = ODE_RungeKutta2(f, x0, result, h, x);
                double exact = Math.Exp(Math.Sin(x));
                if (i % 5 == 0)
                    Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
                x0 = x;
            }
        }
 
        static void TestRungeKutta4()
        {
            double h = 0.001;
            double x0 = 0.0;
            double y0 = 1.0;
            Console.WriteLine("\n Results from the 4th-order Runge-Kutta method with h = {0}\n", h);
            double result = y0;
            for (int i = 0; i < 11; i++)
            {
                double x = 0.1 * i;
                result = ODE_RungeKutta4(f, x0, result, h, x);
                double exact = Math.Exp(Math.Sin(x));
                if (i % 5 == 0)
                    Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
                x0 = x;
            }
        }
 
        static void TestRungeKuttaFehlberg()
        {
            double h = 0.2;
            double x0 = 0.0;
            double y0 = 1.0;
            Console.WriteLine("\n Results from the fourth-order Runge-Kutta-Fehlberg method with h = {0}\n", h);
            double result = y0;
            for (int i = 0; i < 11; i++)
            {
                double x = 0.1 * i;
                result = ODE_RungeKuttaFehlberg(f, x0, result, x, h, 1e-8);
                double exact = Math.Exp(Math.Sin(x));
                if (i % 5 == 0)
                    Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
                x0 = x;
            }
        }
 
        #endregion
 
        static void Main(string[] args)
        {
            TestRungeKutta2();
            TestRungeKutta4();
            TestRungeKuttaFehlberg();
            Console.ReadLine();
        }
    }
}