void NumericJacob (double [], int, double [MAX_DIM][MAX_DIM], double []);

// numerically estimates the jacobian matrix by central finite differences, 
// of the vector forcing function computed by calling CalcFunct 

void NumericJacob (double x_ [], int dim, double J [MAX_DIM][MAX_DIM], double F [])
{
  double work_x [MAX_DIM];
  double F_left [MAX_DIM];
  double F_right [MAX_DIM];
  double delta_x;
  int i, j;
 
 // tunable parameters
  
  double delta_x_rel = 1.0e-7;
  double delta_x_abs = 1.0e-10;


  CalcFunct (x_, dim, F);
  
  for (i = 0; i < dim; i++)
    work_x [i] = x_ [i];
    
  for (i = 0; i < dim; i++)
  {
  	delta_x = delta_x_rel * fabs (x_ [i]);
  	
  	if (delta_x < delta_x_abs) 
      delta_x = delta_x_abs;

    work_x [i] = x_ [i] - delta_x;
    CalcFunct (work_x, dim, F_left);
    
    work_x [i] = x_ [i] + delta_x;
    CalcFunct (work_x, dim, F_right);
    
    work_x [i] = x_ [i];

    for (j = 0; j < dim; j++)
      J [j][i] = (F_right [j] - F_left [j]) / (2 * delta_x);
  }
}