void     
runge_kutta (double *x, double t, double tau,  double *xout){
double F1[NDIM], F2[NDIM], F3[NDIM], F4[NDIM], xtemp[NDIM];
double t_half, half_tau, t_full;
int  i;

  half_tau = 0.5 * tau;
  model (x, t, F1);
  t_half = t + half_tau;
  for (i = 0;i < NDIM; i++)
    xtemp[i] = x[i] + half_tau*F1[i];
  model (xtemp, t_half, F2);
  for (i = 0;i < NDIM; i++)
    xtemp[i] = x[i] + half_tau*F2[i];
  model (xtemp, t_half, F3);
  t_full = t + tau;
  for (i = 0;i < NDIM; i++)
    xtemp[i] = x[i] + tau * F3[i];
  model (xtemp, t_full, F4);
  for (i = 0;i < NDIM; i++)
    xout[i] = x[i] + tau/6.0 *(F1[i] + F4[i] + 2.0*(F2[i] + F3[i]));
}