% Example: three-point approx. of f''; the points are x0-h,x0,x0+h.
f=@(x) exp(x);  x0=1;
% Suppose f values are computed with six decimal digits: 
fa=@(x) round(1e6*f(x))/1e6; % approximate computation of f(x)
df=0.5e-6; % estimation of the error in f
ddf_e=exp(x0); % exact f''
% What is the best h? Experiments with different h:
hh=2.^(4:-0.25:-18); 
N=length(hh);err=zeros(1,N);
for i=1:N
    h=hh(i);
    df_a=(fa(x0-h)-2*fa(x0)+fa(x0+h))/h^2; % approximate f''
    err(i)=abs(df_a-ddf_e); % error in f''
end
% error estimate: 
k=exp(1)/12; h2=hh.^2; err_est=4*df./h2+k*h2;
% results:
loglog(hh,err,'b',hh,err_est,'--r','LineWidth',2);
legend('real error','estimate')
grid on;xlabel('h');ylabel('error');axis tight
title('Example: numerical calculation of d^2f/dx^2')
[e_min,i]=min(err);
he=(48*df/exp(1))^0.25; ee=4*df/he^2+exp(1)*he^2/12;
fprintf('Experiment: the best h = %4.2e   min error =  %4.2e\n',hh(i),e_min);
fprintf('Estimate:   the best h = %4.2e   min error =  %4.2e\n',he,ee);