function [h,ef,ec,ei] = fordifferr( f, a, dfa, k )
  function f = fibonacci(n)
    f = zeros(1,n);
    f(1) = 1;
    f(2) = 1;
    for i=3:n,
        f(i) = f(i-1)+f(i-2);
    end;
  end
  ns = fibonacci(k+2);
  h  = 1.0./ns(3:end); % Throw away 1's at start
  x  = a + h;
  ef = zeros(1,k);
  ec = zeros(1,k);
  ya = feval( f, a );
  ef = ( feval( f, x ) - ya )./h/dfa - 1 ;
  ec = ( feval( f, a+h/2 ) - feval( f, a-h/2 ) )./h/dfa - 1;
  ei = ( imag( feval( f, a+i*h ) )./h/dfa ) - 1;
end
   