function y = dx_central(x,h,f)

  y = ( feval(f,x+h) - feval(f,x-h) )/(2*h);