function num_dif3
%

f01 = @(x) x.^9;
f01d= @(x) 9*x.^8;

x = 1; 
N = 17;
h = 1; 
hvals = zeros(1,N);
dfbydx = hvals;

exact = f01d(x); 

disp(' eval(diff(x^9, x), x=1.0) = 9 ' )
disp('   h    Exact  Computed   Error')

for i = 1:N
  h = h/10; 
  hvals(i) = h; 
  dfbydx(i) = (f01(x+h) - f01(x))/h;

  fprintf('%6.1e   9 %11.2e %10.2e \n', ...
          hvals(i), dfbydx(i), exact-dfbydx(i))
end;

loglog(hvals, abs(dfbydx - exact)/exact, '*', 'MarkerSize', 8);
axis([1e-18 1 1e-8 5])
xlabel('h value'); 
ylabel('error in approximation'); 
grid

end % function
