C#龙格库塔解微分方程

请问大家，C#解龙格库塔该怎么写，刚学一点C#基本语句，只会定义k的四行，不知道完整的该怎么写

float dydx(float x, float y)
{
    return ((x - y) / 2);
}
float rungeKutta(float x0, float y0, float x, float h)
{
    // Count number of iterations using step size or
    // step height h
    int n = (int)((x - x0) / h);

    float k1, k2, k3, k4, k5;

    // Iterate for number of iterations
    float y = y0;
    for (int i = 1; i <= n; i++)
    {
        // Apply Runge Kutta Formulas to find
        // next value of y
        k1 = h * dydx(x0, y);
        k2 = h * dydx(x0 + 0.5f * h, y + 0.5f * k1);
        k3 = h * dydx(x0 + 0.5f * h, y + 0.5f * k2);
        k4 = h * dydx(x0 + h, y + k3);

        // Update next value of y
        y = y + (1.0f / 6.0f) * (k1 + 2 * k2 + 2 * k3 + k4); ;

        // Update next value of x
        x0 = x0 + h;
    }

    return y;
}

//调用：
rungeKutta(1, 2, 3, 4);

/*函数引用自：
https://www.geeksforgeeks.org/runge-kutta-4th-order-method-solve-differential-equation/
*/

题主可以试一下。

using System;
using System.Collections;
using System.Collections.Generic;
 
namespace Legalsoft.Truffer.Algorithm
{
    /// <summary>
    /// 给定微分方程的一阶偏导方程
    /// </summary>
    /// <param name="x"></param>
    /// <param name="y"></param>
    /// <returns></returns>
    public delegate double SDE_Equation(double x, double y);
 
    /// <summary>
    /// 解微分方程的龙格-库塔四阶方法
    /// </summary>
    public static partial class Algorithm_Gallery
    {
        public static SDE_Equation dydx = null;
 
        /// <summary>
        /// Finds value of y for a given x
        /// using step size h
        /// and initial value y0 at x0.
        /// </summary>
        /// <param name="x0">初值</param>
        /// <param name="y0">初值</param>
        /// <param name="x">求值点</param>
        /// <param name="h">步长</param>
        /// <returns></returns>
        public static double Runge_Kutta_4th_Order(double x0, double y0, double x, double h)
        {
            int n = (int)((x - x0) / h);
 
            double y = y0;
            for (int i = 1; i <= n; i++)
            {
                double k1 = h * (dydx(x0, y));
                double k2 = h * (dydx(x0 + 0.5 * h, y + 0.5 * k1));
                double k3 = h * (dydx(x0 + 0.5 * h, y + 0.5 * k2));
                double k4 = h * (dydx(x0 + h, y + k3));
 
                y = y + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
 
                x0 = x0 + h;
            }
 
            return y;
        }
    }
}