function [d] = richardson(f, x, n, h0)
% richardson extrapolation for derivatives
% f is the function handle of the function
% x is the value of x for f'(x)
% n is number of rows in table
% h0 is the initial step size

% This is a conversion of the Cheney and Kincaid
% 0-based algorithm to a 1-based algorithm

d = zeros(n,n);
h = h0;
for i = 0 : n-1 % rows 1 to n of the table
 
   d(i+1,1) = (f(x + h) - f(x - h))/(2*h); % column 1
   
   % Use richardson extrapolation across row i+1
   % to the diagonal
   
   for j = 1:i
      d(i+1,j+1) = d(i+1,j) + (d(i+1,j) - d(i,j)) / (4^j-1);
   end
   h = h/2;
end