function [d,err] = Derivative2(fname,a,delta,M2,M3)
% d = Derivative(fname,a,delta,M2,M3);
% fname  a string that names a function f(x) whose derivative
% at x = a is sought. delta is the absolute error associated with
% an f-evaluation and M2 is an estimate of the second derivative
% magnitude near a, and M3 is an estimate of the third derivative.
% d is an approximation to f'(a) using the better of forward or
% central divided difference, and err is an estimate 
% of its absolute error.
%
% Usage:
%     [d,err] =  Derivative(fname,a)   
%     [d,err] =  Derivative(fname,a,delta) 
%     [d,err] =  Derivative(fname,a,delta,M2)
%     [d,err] =  Derivative(fname,a,delta,M2,M3)
%
% Caveat:  The original function Derivative didn't do any
% error checking (e.g. delta, M2, and M3 should be positive),
% so neither are we.

if nargin <= 4
   % for no good reason, guess that
   % third derivative bound is 1.
   M3 = 1;
end
if nargin <= 3
   % No derivative bound supplied, so assume the
   % second derivative bound is 1.
   M2 = 1;
end
if nargin == 2
   % No function evaluation error supplied, so
   % set delta to eps.
   delta = eps;
end
% Which method is better?
if 81*delta*M3^2 < 256*M2^3

% Compute optimum h and central divided difference
  hopt = (6*delta/M3)^(1/3);
  d = (feval(fname,a+hopt) - feval(fname,a-hopt))/2/hopt;
  err = 3*(M3*delta^2/6)^(1/3);

else

% Compute optimum h and forward divided difference
  hopt = 2*sqrt(delta/M2);
  d   = (feval(fname,a+hopt) - feval(fname,a))/hopt;
  err = 2*sqrt(delta*M2);

end