Below is the program illustrating R-K Method in C#.

Code:
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Text;
using System.Windows.Forms;

namespace RKMethod
{
    public partial class Form1 : Form
    {
        public enum Variables
        {
            X, Y,CONSTANT
        }
        struct Polynomial
        {
            public float Exponent;
            public float Coefficient;
            public Variables Variable;
        };

        private Polynomial[] Poly;
        private int polyCount = 0;

        public Form1()
        {
            InitializeComponent();

            polyCount = 0;

            initPoly(2);

            /*This polynomial is equal to
             * x^2 + Y
             */
 
            addTerm(Variables.X, 3, 1);
            addTerm(Variables.Y, 0.5f, 1);
        }

        public void initPoly(int termCount) {

            Poly = new Polynomial[termCount];
        }
        public void addTerm(Variables v, float Coefficient, float Exponent)
        {
            Poly[polyCount].Coefficient = Coefficient;
            Poly[polyCount].Exponent = Exponent;
            Poly[polyCount].Variable = v;

            polyCount++;
        }
        public float Fxy(float xValue,float yValue)
        {
            float val = 0;
            float e = 1;

            for (int i = 0; i < polyCount; i++)
            {
                bool isVar = false;
                float value;

                if (Poly[i].Variable == Variables.X)
                {
                    isVar = true;
                    value = xValue;
                }
                else if (Poly[i].Variable == Variables.Y)
                {
                    isVar = true;
                    value = yValue;
                }
                else
                {
                    isVar = false;
                    value = -1;
                }

                if (isVar)
                {
                    for (int j = 0; j < Poly[i].Exponent; j++)
                        e *= value;

                    e *= Poly[i].Coefficient;
                    val += e;
                }
                else
                {
                    val += (float)Poly[i].Coefficient;
                }
                e = 1;
            }   
            return val;
        }
        public float RK(float X,float Y,float newX) {

            float k1, k2, k3, k4, k, h, newY;

            h = newX - X;
            k1 = h * Fxy(X, Y);
            k2 = h * Fxy(X + (h / 2), Y + (k1 / 2));
            k3 = h * Fxy(X + (h / 2), Y + (k2 / 2));
            k4 = h * Fxy(X + h, Y + k3);

            k = ((k1 + (2 * k2) + (2 * k3) + k4) / 6);
            
            newY = Y + k;

            return newY;
        }
        private void button1_Click(object sender, EventArgs e)
        {
            float v = RK(0, 1, 0.1f);
            MessageBox.Show(v.ToString());
        }
    }
}