Numerical implementation of n-th derivative of f(x)? [closed]

I implemented a C++ code to numerically solve the n-th derivative of a function in a point x_0:

double n_derivative( double ( *f )( int, double ), double x_0, int n )
 {
  if( n == 0 ) return f( x_0 );
  else 
   {
    const double h = pow( __DBL_EPSILON__, 1/3 );
    double x_1 = x_0 - h;
    double x_2 = x_0 + h;
    double first_term = n_derivative( f, x_2, n - 1 );
    double second_term = n_derivative( f, x_1, n - 1);

    return ( first_term - second_term ) / ( 2*h );
   }
 }

I was wondering if this is for you a good implementation or if there can be a way to better write it in C++. The problem is that I noticed that the n-th derivative diverges for values of n higher than 3. Do you know how to solve this?

It is not a good implementation

At least these problems.

Integer math

Use FP math as 1/3 is zero.

1/3 --> 1.0/3

Using the cube root optimal for n==1

But not certainly other n. @Eugene

Wrong epsilon

Below code is only useful for |x_0| about 1.0. When x_0 is large, x_0 - h may equal x_0. When x_0 is small, x_0 - h may equal -h.

OP's +/- some epsilon is good for fixed point, but double is a floating point.

// Bad
const double h = pow( __DBL_EPSILON__, 1.0/3 );
double x_1 = x_0 - h;

A relative scaling is needed.

#define EPS cbrt(DBL_EPSILON) // TBD code to well select this 
if (fabs(x_0) >= DBL_MIN && isfinite(x_0)) {
  double x_1 = x_0*(1.0 - EP3);
  double x_2 = x_0*(1.0 + EPS);
  double h2 = x_2 - x_1;
  ...
} else {
  TBD_Code for special cases
}

Invalid code

f is double ( *f )( int, double ), but call is f( x_0 )

Minor: confusing names

Why first_term with x_2 and second_term with x_1?