clear;
close all;
format long;
N=10;
x=1;
for ii=1:N
h_loop(ii)=2^(-ii);
end

% table headings:
disp(' ')
disp('       h              D1             D2')

for ii=1:length(h_loop)
    h=h_loop(ii);

D1u=(fun(x+h)-2*fun(x)+fun(x-h))/(h^2);
E1u(ii)=D1u-d2fun(1);
%Dmu=(fun(x)-fun(x-h))/h;
%Emu(ii)=Dmu-dfun(1);
D2u=(fun(x+2*h)-2*fun(x)+fun(x-2*h))/((2*h)^2);
E2u(ii)=D2u-d2fun(1);

% print line of table:
   disp(sprintf(' %13.4e  %13.6e  %13.4e  ',...
       h,E1u(ii),E2u(ii)));
end
loglog(h_loop,abs(E1u),'.-');
hold on;
loglog(h_loop,abs(E2u),'o-');

       legend('D1','D2');
       xlabel('h')
       ylabel('Error')
       title('u^{2}(x)')