[D,err,relerr,n]=diffext(f,x,delta,toler)

函数 diffext 用于数值求解函数 f 在点 x 处的一阶导数。它使用了外推型 Richardson 外推公式，并自适应调整步长以达到所需的精度。

参数 delta 是初始步长，toler是允许的误差限度。函数返回三个值：估计的导数值 D、误差 Err 和相对误差 RelErr，以及使用的步长数 N。

以下是函数 diffext 的 MATLAB 代码实现：

function [D,Err,RelErr,N] = diffext(f, x, delta, toler)
% DIFFEXT Numerical differentiation using Richardson extrapolation.
%   D = DIFFEXT(F, X, DELTA) estimates the first derivative of F at X
%   using a central difference formula with step size DELTA.
%   The function F should accept a scalar input and return a scalar output.
%
%   D = DIFFEXT(F, X, DELTA, TOLER) uses the given tolerance TOLER
%   to control the accuracy of the estimate. The algorithm adapts
%   the step size so that the relative error is less than TOLER.
%
%   [D, ERR, RELERR, N] = DIFFEXT(...) also returns the estimated error
%   ERR, the relative error RELERR, and the number of steps N used.
%
%   Reference: W. H. Press, et al., "Numerical Recipes in C," 2nd ed., Sec. 5.7.

% Initializations
D1 = (f(x+delta)-f(x-delta))/(2*delta);
N = 1;      % Number of steps
delta = delta/2;
D2 = (f(x+delta)-f(x-delta))/(2*delta);
N = N + 1;
Err = abs(D2 - D1)/3;
RelErr = Err/max(abs(D1), abs(D2));

% Adapt step size
while RelErr > toler
    D1 = D2;
    delta = delta/2;
    D2 = (f(x+delta)-f(x-delta))/(2*delta);
    N = N + 1;
    Err = abs(D2 - D1)/3;
    RelErr = Err/max(abs(D1), abs(D2));
end

% Final estimate of derivative
D = D2;

例子：

% Define the function f(x) = x^3
f = @(x) x.^3;

% Estimate the derivative of f at x = 1 using step size delta = 0.1
x = 1;
delta = 0.1;
[D, Err, RelErr, N] = diffext(f, x, delta);

% Display the results
fprintf('The estimated derivative of f at x = %g is %g\n', x, D);
fprintf('The error estimate is %g\n', Err);
fprintf('The relative error is %g\n', RelErr);
fprintf('The number of steps used is %d\n', N);