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

% table headings:
disp(' ')
disp('       h              D+             D-             D0')

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

Dpu=(fun(x+h)-fun(x))/h;
Epu(ii)=Dpu-dfun(1);
Dmu=(fun(x)-fun(x-h))/h;
Emu(ii)=Dmu-dfun(1);
Dou=(fun(x+h)-fun(x-h))/(2*h);
Eou(ii)=Dou-dfun(1);
% print line of table:
 disp(sprintf(' %13.4e  %13.6e  %13.4e  %13.4e ',...
       h,Epu(ii),Emu(ii),Eou(ii)));
end

%plot(log10(h_loop),log10(abs(Epu)),'.-');
loglog(h_loop,abs(Epu),'.-');
hold on;
%plot(log10(h_loop),log10(abs(Eou)),'o-');
loglog(h_loop,abs(Eou),'o-');

       legend('D+','D0');
       xlabel('h')
       ylabel('Error')
       title('example 1')