f = inline('exp(x)', 'x');
 fprime = inline('exp(x)','x');

 x0 = 1;
 h = 2.^-[0:20];               % h = 1, 1/2, 1/4, ...
 phi = zeros(length(h),1);
 psi = zeros(length(h)-1,1);
 theta = zeros(length(h)-2,1);
 exact = fprime(x0);

 fprintf('       h               phi(h)             phi(h)           theta(h)       error in theta \n')
 fprintf('  ------------------------------------------------------------------------------------------\n')
 for j=1:length(h)
    phi(j) = (f(x0+h(j)) - f(x0))/h(j);
 end
 for j=1:length(h)-1
    psi(j) = 2*phi(j+1)-phi(j);
 end

 for j=1:length(h)-2
    theta(j) = (4*psi(j+1)-psi(j))/3;
    fprintf('   %10.5e    %14.12f     %14.12f     %14.12f    %10.5e \n', ...
             h(j), phi(j), psi(j), theta(j), abs(theta(j)-exact))
 end 

 pause
 figure(1), clf
 phiplt = loglog(h, abs(phi-exact), 'b-', 'linewidth',2); hold on
          loglog(h, abs(phi-exact), 'b.', 'markersize', 20)
 psiplt = loglog(h(1:end-1), abs(psi-exact), 'k-', 'linewidth',2);
          loglog(h(1:end-1), abs(psi-exact), 'k.', 'markersize', 20)
 theplt = loglog(h(1:end-2), abs(theta-exact), 'm-', 'linewidth',2);
          loglog(h(1:end-2), abs(theta-exact), 'm.', 'markersize', 20)
 set(gca,'fontsize',20)
 xlabel('h', 'fontsize', 20)
 ylabel('absolute error', 'fontsize', 20)
 axis([1e-7 1e1 1e-15 1e2]) 
 set(gca,'XDir','reverse')
 legend([phiplt psiplt theplt], '\phi(h)', '\psi(h)', '\theta(h)')
 pause
 loglog(h, h, 'r-', 'linewidth', 2)
 loglog(h, h.^2, 'r-', 'linewidth', 2)
 loglog(h, h.^3, 'r-', 'linewidth', 2)

 pause

 format long
 R = zeros(3);
 R(:,1) = phi(1:3);
 R(2:3,2) = psi(2:3);
 R(3,3)  = theta(3);
 fprintf('\n\n\n')
 R