function d = richardson( f, x, n, h )
  % ----------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 9, 1999
  %
  % Usage: D = richardson(@f, x, n, h)
  %
  % Returns:
  %	two-dimensional array of extrapolation values.  The n+1,n+1 value
  %	of this array should be the most accurate.
  %
  % Parameters:
  %	f:	handle of function to find derivative of
  %	x:	value of x to find derivative at
  %	n:	number of levels of extrapolation
  %	h:	initial stepsize
  %
  % This function implements Richardson's Extrapolation to determine an
  % approximation of f'(x) at a particular value of x.
  %
  % Based on an algorithm in "Numerical Mathematics and Computing"
  % 4th Edition, by Cheney and Kincaid, Brooks-Cole, 1999.
  %
  % Note: There are several changes that could be made to make this code
  % more efficient.  However, it follows the algorithm in the text very
  % closely.  The main difference is that arrays in octave start with a
  % subscript of 1 rather than 0.
  % ----------------------------------------------------------------------

  % d(n+1,n+1) will contain the most accurate approximation to f'(x).
  d = zeros( n + 1, n + 1 );
  for i = 1 : n + 1
    d(i,1) = 0.5 * ( f( x + h ) - f( x - h ) ) / h;
    pow4 = 4;
    for j = 1 : i - 1
      d(i,j+1) = d(i,j) + ( d(i,j) - d(i-1,j) ) / ( pow4 - 1 );
      pow4 = 4 * pow4;
    end
    h = 0.5 * h;
  end
end
