Solving a differential equation

abstract public class FirstDerivativeBase
{
    abstract public double GetValue(double t, double y);
}

class CircuitV : FirstDerivativeBase
{
    double c;
    double r;

    public CircuitV(double c, double r)
    {
       this.c = 1;
       this.r = 1;

       if (c > 0)
          this.c = c;
       if (r > 0)
          this.r = r;
    }

    public override double GetValue(double t, double q)
    {
       double u = q / c;
       double i = -u / r;
       return i;
    }
}

public class IntegratorRK4V
{
    FirstDerivativeBase derivative;
    double t0, y0, h;
    double d0;

    public IntegratorRK4V(FirstDerivativeBase derivative,
        double t0, double y0, double h)
    {
       this.derivative = derivative;
       this.t0 = t0;
       this.y0 = y0;
       this.h = h;
       this.d0 = derivative.GetValue(t0, y0);
    }

    public int Step(out double t1, out double y1, out double d1)
    {
       double k1 = d0;
       double k2 = derivative.GetValue(t0 + h / 2, y0 + h / 2 * k1);
       double k3 = derivative.GetValue(t0 + h / 2, y0 + h / 2 * k2);
       double k4 = derivative.GetValue(t0 + h, y0 + h * k3);

       y1 = y0 + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
       t1 = t0 + h;
       d1 = derivative.GetValue(t1, y1);

       t0 = t1;
       y0 = y1;
       d0 = d1;

       return 0;
    }
}

static void Test1()
{
    CircuitV circV = new CircuitV(0.1,500);

    double t0 = 0;
    double tmax = 90;
   double q0 = 25;
    double h = 5;
    double i0 = circV.GetValue(t0, q0);
    IntegratorRK4V integratorRK4V = new IntegratorRK4V(circV,t0,q0,h);

    Console.WriteLine("TEST1, virtual function");
    Console.WriteLine("t,s     Q,C     I,A");
    Console.WriteLine("-------------------");
    Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}",t0,q0,-i0);

    double t, q, i;
    do
    {
       integratorRK4V.Step(out t, out q, out i);
       Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t, q, -i);
    } while (t < tmax);
    Console.WriteLine();
}

class IntegratorRK4I
{
    FirstDerivative firstDerivative;
    double t0, y0, h;
    double d0;

    public IntegratorRK4I(
        FirstDerivative firstDerivative, double t0, double y0,
        double h, double d0)
    {
       this.t0 = t0;
       this.y0 = y0;
       this.h = h;
       this.d0 = d0;
       this.firstDerivative = firstDerivative;
    }

    public int Step(out double t1, out double y1, out double d1)
    {
       double k1 = d0;
       double k2 = firstDerivative.GetValue(t0 + h / 2, y0 + h / 2 * k1);
       double k3 = firstDerivative.GetValue(t0 + h / 2, y0 + h / 2 * k2);
       double k4 = firstDerivative.GetValue(t0 + h, y0 + h * k3);

       y1 = y0 + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
       t1 = t0 + h;
       d1 = firstDerivative.GetValue(t1, y1);

       t0 = t1;
       y0 = y1;
       d0 = d1;

       return 0;
    }
}

static void Test2()
{
    CircuitI circI = new CircuitI(0.1, 500);
    FirstDerivative firstDerivative = (FirstDerivative)circI;

    double t0 = 0;
    double tmax = 90;
    double q0 = 25;
    double h = 5;
    double i0 = firstDerivative.GetValue(t0, q0);
    IntegratorRK4I integratorRK4I = new IntegratorRK4I(
        firstDerivative, t0, q0, h, i0);

    Console.WriteLine("TEST2, interface");
    Console.WriteLine("t,s     Q,C     I,A");
    Console.WriteLine("-------------------");
    Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t0, q0, -i0);

    double t, q, i;
    do
    {
       integratorRK4I.Step(out t, out q, out i);
       Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t, q, -i);
    } while (t < tmax);
    Console.WriteLine();
}

class IntegratorRK4D
{
    public delegate double GetDerivative(double t, double y);
    GetDerivative getDerivative;
    double t0, y0, h;
    double d0;

    public IntegratorRK4D(GetDerivative getDerivative, double t0,
        double y0, double h, double d0)
    {
       this.getDerivative = getDerivative;
       this.t0 = t0;
       this.y0 = y0;
       this.h = h;
       this.d0 = d0;
    }

    public int Step(out double t1, out double y1, out double d1)
    {
       double k1 = d0;
       double k2 = getDerivative(t0 + h / 2, y0 + h / 2 * k1);
       double k3 = getDerivative(t0 + h / 2, y0 + h / 2 * k2);
       double k4 = getDerivative(t0 + h, y0 + h * k3);

       y1 = y0 + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
       t1 = t0 + h;
       d1 = getDerivative(t1, y1);

       t0 = t1;
       y0 = y1;
       d0 = d1;

       return 0;
    }
}

static void Test3()
{
    CircuitD circD = new CircuitD(0.1, 500);

    double t0 = 0;
    double tmax = 90;
    double q0 = 25;
    double h = 5;
    double i0 = circD.GetDerivative(t0, q0);
    IntegratorRK4D integratorRK4D = new IntegratorRK4D(circD.GetDerivative,
        t0, q0, h, i0);

    Console.WriteLine("TEST3, delegate");
    Console.WriteLine("t,s     Q,C     I,A");
    Console.WriteLine("-------------------");
    Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t0, q0, -i0);

    double t, q, i;
    do
    {
       integratorRK4D.Step(out t, out q, out i);
       Console.WriteLine("{0,3:F0}{1,8:F1}{2,8:F2}", t, q, -i);
    } while (t < tmax);
    Console.WriteLine();
}

class FlyingBall : SecondDerivative
{
    double dB;    // Ball diameter, m
    double roB;   // Ball density, kg/m^3
    double cD;    // Drag coefficient
    double g;     // Gravitational acceleration, m/s^2
    double roA;   // Air density
    double mB;    // Ball mass, kg
    double aB;    // Ball cross-section area, m^2

    public FlyingBall()
    {
       dB = 0.1;
       roB = 600;
       cD = 0.1;
       g = 9.81;
       roA = 1.29;

       double v = Math.PI * Math.Pow(dB, 3) / 6;
       mB = roB * v;
       aB = 0.25 * Math.PI * dB * dB;
    }

    double SecondDerivative.GetValue(double t, double y, double v)
    {
       double f = -mB * g -Math.Sign(v) * 0.5 * cD * roA * v * v * aB;
       double d2y = f / mB;
       return d2y;
    }
}

class IntegratorRKN
{
    SecondDerivative secondDerivative;
    double t0, y0, dy0, h;
    double d2y0;

    public IntegratorRKN(SecondDerivative secondDerivative, double t0,
        double y0, double dy0, double h, double d2y0)
    {
       this.secondDerivative = secondDerivative;
       this.t0 = t0;
       this.y0 = y0;
       this.dy0 = dy0;
       this.h = h;
       this.d2y0 = d2y0;
    }

    public int Step(out double t1, out double y1, out double dy1,
        out double d2y1)
    {
       double h2 = h * h;

       double k1 = d2y0;
       double k2 = secondDerivative.GetValue(
           t0 + h/2, y0 + h/2 * dy0 + h2/8 * k1, dy0 + h/2 * k1);
       double k3 = secondDerivative.GetValue(
           t0 + h/2, y0 + h/2 * dy0 + h2/8 * k2, dy0 + h/2 * k2);
       double k4 = secondDerivative.GetValue(
           t0 + h,   y0 +   h * dy0 + h2/2 * k3, dy0 +   h * k3);

       t1 = t0 + h;
       y1 = y0 + h * dy0 + h2/6 * (k1 + k2 + k3);
       dy1 = dy0 + h/6 * (k1 + 2 * k2 + 2 * k3 + k4);
       d2y1 = secondDerivative.GetValue(t1, y1, dy1);

       t0 = t1;
       y0 = y1;
       dy0 = dy1;
       d2y0 = d2y1;

       return 0;
}

static void Test4()
{
    FlyingBall ball = new FlyingBall();
    SecondDerivative acceleration = (SecondDerivative) ball;

    double t0 = 0;
    double tmax = 10;
    double y0 = 0;
    double h = 0.5;
    double v0 = 50;
    double a0 = acceleration.GetValue(t0, y0, v0);
    IntegratorRKN integratorRKN = new IntegratorRKN(acceleration, t0, y0, v0,
        h, a0);

    Console.WriteLine("TEST4, equation of motion");
    Console.WriteLine(" t,s     y,m   v,m/s  a,m/s2");
    Console.WriteLine("----------------------------");
    Console.WriteLine("{0,4:F1}{1,8:F2}{2,8:F2}{3,8:F2}", t0, y0, v0, a0);

    double t, y, v, a;
    do
    {
       integratorRKN.Step(out t, out y, out v, out a);
       Console.WriteLine("{0,4:F1}{1,8:F2}{2,8:F2}{3,8:F2}", t, y, v, a);
    } while (t < tmax);
    Console.WriteLine();
}