Интегрирование дифференциальных уравнений движения для орбиты в Рунге Кутте в C

Question

Я пытаюсь интегрировать дифференциальные уравнения для планетарной орбиты, используя runge-kutta четвертого порядка в C. Первое уравнение, которое я попытался интегрировать, было для позиции, dr/dt = sqrt((2/mu)*(E-(k/r)) - pow(1, 2)/(pow(mu, 2)*pow(r, 2))). Программа компилируется правильно, но продолжает возвращать nan, я довольно плохо знаком с C, поэтому мне интересно, почему это может быть. Исходный код, который я использую, является следующим:

    double rk4 ( double t0, double u0, double dt, double f ( double t, double u ) ) 
    {


    /* Parameters:

    Input, double T0, the current time.

    Input, double U0, the solution estimate at the current time.

    Input, double DT, the time step.

    Input, double F ( double T, double U ), a function which evaluates
    the derivative, or right hand side of the problem.

    Output, double RK4, the fourth-order Runge-Kutta solution estimate
    at time T0+DT.
*/
    double f0;
    double f1;
    double f2;
    double f3;
    double t1;
    double t2;
    double t3;
    double u;
    double u1;
    double u2;
    double u3;
     /*
  Get four sample values of the derivative.
*/
  f0 = f ( t0, u0 );

  t1 = t0 + dt / 2.0;
  u1 = u0 + dt * f0 / 2.0;
  f1 = f ( t1, u1 );

  t2 = t0 + dt / 2.0;
  u2 = u0 + dt * f1 / 2.0;
  f2 = f ( t2, u2 );

  t3 = t0 + dt;
  u3 = u0 + dt * f2;
  f3 = f ( t3, u3 );
/*
  Combine them to estimate the solution.
*/
  u = u0 + dt * ( f0 + 2.0 * f1 + 2.0 * f2 + f3 ) / 6.0;

  return u;
}

И мой код выглядит следующим образом:

    #include <stdio.h>
    #include <math.h>
    #include "rk4.c"

    double f ( double t, double u );
    double rk4 ( double t0, double u0, double dt, double f ( double t, double u ));

    int main(void)
    {
    double t0 = 0, u0 = 1.0167103, dt = 0.01, E, l = 2.663e40, k = -1.32754125e20, mu = 5.9719821e24, G = 6.67408e-11, M = 1.989e30, h, e, u, t; 

    printf("%lf\n", rk4 ( t0, u0, dt, f));

    return 0;
    }

    double f ( double t, double u ) 
    {
    double E, l = 2.663e40, k = -1.32754125e20, u0 = 1.0167103, mu = 5.9719821e24, e = 0.0167, h;

    h = l/mu;
    E = ((pow(e, 2)*pow(mu, 3) - pow(mu, 3))/2*pow(h, 2));

    return sqrt((2/mu)*(E-(k/u)) - pow(l, 2)/((pow(mu, 2))*(pow(u, 2))));
    }