© 2021 Neveles. All rights reserved.

int funForRKF(int n, double t, double* x, double* dxdt)
{
  dxdt[0] = x[1]; // m1' = m2
  dxdt[2] = x[3]; // m3' = m4

  dxdt[1] = -(1 - u) * x[0] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
    (x[0] - 1) / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));

  dxdt[3] = -(1 - u) * x[2] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
    x[2] / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));

  return 0;
}