  function d = Derivative(fname,a,delta,M2);
%
% Pre:
%    fname  a string that names an available function.
%    a:     the x-value at which an estimate of f'(x) is required.
%    delta: absolute error in function evaluation.
%    M2:    an estimate of the second derivative magnitude near a.
%
% Post:
%    d:    an approximation of f'(a).
%    err:  an estimate of the error in d.
%
% Usage:
%     d =  Derivative(fname,a)   
%     d =  Derivative(fname,a,delta) 
%     d =  Derivative(fname,a,delta,M2)

   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

% Compute optimum h and divided difference

   hopt = 2*sqrt(delta/M2);
   d   = (feval(fname,a+hopt) - feval(fname,a))/hopt;

      
      
