function [D,err,relerr,n] = diffext(f,x,delta,toler)
% [D,err,relerr,n] = diffext(f,x,delta,toler)
% Numerical approximation for f'(x).
% The method is Richardson`s extrapolation.
% f is the function, input.
% x is differentiation point, input.
% delta is the error goal, input.
% toler is the relative error goal, input.
% D is the matrix of approximate derivatives, output.
% error is the error bound, output.
% relerr is the relative error bound, output.
% n is the coordinate of the "best approximation", output.
err = 1;
relerr = 1;
h = 1;
j = 1;
D(0+1,0+1) = (feval(f,x+h) - feval(f,x-h))/(2*h);
while relerr>toler & err>delta & j<12
h = h/2;
D(j+1,0+1) = (feval(f,x+h) - feval(f,x-h))/(2*h);
for k = 1:j,
D(j+1,k+1) = D(j+1,k-1+1) + (D(j+1,k-1+1)-D(j-1+1,k-1+1))/(4^k -1);
end
err = abs(D(j+1,j+1)-D(j-1+1,j-1+1));
relerr = 2*err/(abs(D(j+1,j+1))+abs(D(j-1+1,j-1+1))+eps);
j = j+1;
end
[n n] = size(D);