// Runge-Kutta - 4th order
    void RungeKutta()
    {
        // The stepsize
        const double h = 0.005;
        
        Vecd k1v = h * gVelocities;
        Vecd k1f = h * SumForces(gPositions, gVelocities);
        
        Vecd k2v = h * (gVelocities + k1f/2.0);
        Vecd k2f = h * SumForces(gPositions + k1v/2.0, gVelocities + k1f/2.0);
        
        Vecd k3v = h * (gVelocities + k2f/2.0);
        Vecd k3f = h * SumForces(gPositions + k2v/2.0, gVelocities + k2f/2.0);
        
        Vecd k4v = h * (gVelocities + k3f);
        Vecd k4f = h * SumForces(gPositions + k3v, gVelocities + k3f);


        gVelocities += k1f/6 + k2f/3 + k3f/3 + k4f/6;
        gPositions  += k1v/6 + k2v/3 + k3v/3 + k4v/6;
    }