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

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) = (f401(x+h) - f401(x))/h;
  
  %disp(sprintf('%6.1e   9 %11.2e %10.2e', hvals(i), dfbydx(i), 9-dfbydx(i)))
  fprintf('%6.1e   9 %11.2e %10.2e \n', hvals(i), dfbydx(i), 9-dfbydx(i))
end;

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


function f = f401(x)
f = x.^9;

function f = f401_d(x)
f = 9*x.^8;
