function y = CentDiff( f,x,h )
% CentDiff(f,x,h) returns the central difference
%   approximation for the derivative of the
%   function f at x, using step size h > 0,
%   with the formula (f(x+h/2) - f(x-h/2)/h).
  y = 0;
  if(h <= 0)
    fprintf('h must be positive');
  else
    y = feval(f,x+h/2) - feval(f,x-h/2);
    y = y/h;
  end
end