clear
x=1;
h=1.e-1;
disp(' ')
disp('       h              Dpu            Dmu             D0u             D3u')
   
Dp=cos(x);
hvals = [1e-1 5e-2 1e-2 5e-3 1e-3];
Epu = [];  Emu = [];  E0u = [];  E3u = [];
Opu = [];  Omu = [];  O0u = [];  O3u = [];
for i=1:length(hvals)
    h = hvals(i);
%i=1;
Dpu=(sin(x+h)-sin(x))/h;
Dmu=(sin(x)-sin(x-h))/h;
D0u=(sin(x+h)-sin(x-h))/(2*h);
D3u=(2*sin(x+h)+3*sin(x)-6*sin(x-h)+sin(x-2*h))/(6*h);
Epu(i)=Dpu-Dp;
Emu(i)=Dmu-Dp;
E0u(i)=D0u-Dp;
E3u(i)=D3u-Dp;
%end
% table headings:
disp(sprintf('%13.4e   %13.4e   %13.4e   %13.4e  %13.4e',...
                 h,Epu(i),Emu(i),E0u(i),E3u(i)))
end
disp(' ')
disp('       h              Dpu            Order')
Opu(1)=0;
disp(sprintf('%13.4e   %13.4e   %13.3f  ',...
                 hvals(1),Epu(1),Opu(1)))
for i=2:length(hvals)
Opu(i)=log(Epu(i)/Epu(1))/log(hvals(i)/hvals(1));
disp(sprintf('%13.4e   %13.4e   %13.3f  ',...
                 hvals(i),Epu(i),Opu(i)))
end
figure(1);
clf
loglog(hvals,abs(Epu),'o-')
axis([5e-4 .2 1e-12 1]);
hold on
loglog(hvals,abs(E0u),'d-')
loglog(hvals,abs(E3u),'s-')
legend('Dpu','D0u','D3u');
hold off
%loglog(hvals,abs(Emu),'o-')

logh=log10(hvals);
[c,d]=polyfit(logh,log10(Epu),1);
disp(sprintf('Order = %13.3f',c(1)))
figure(2);
clf;
set(axes,'fontsize',18)
hold on
plot(logh,log10(Epu),'s');
% plot the data points
err_fit=c(1)*logh+c(2);
plot(logh,err_fit);
%draw the line that fits the data in a least-square sense
xlabel('-log_{10} h');
ylabel('log_{10} err');
hold off

