Programa 6.2 (Derivacion numerica mediante extrapolacion - pagina 349).

%Programa 6.2 (Derivacion numerica mediante extrapolacion - pagina 349).
% Construccion de la tabla D(j,k) (con k<=j) de aproximaciones numericas a f'(x) en la que se
% utiliza f'(x)=D(n,n) como respuesta final. Las aproximaciones D(j,k) se almacenan en una
% matriz triangular inferior cuya primera columna viene dada por
%
% D(j,k) = ( ( f(x+c) - f(x-c) ) / (2·c) ), con c=2^(-j)·h
% y cuyos elementos en la fila j-esima, para j>=2, son
% D(j,k) = D(j,k-1) + ( D(j,k-1) - D(j-1,k-1) ) / (4^k - 1), donde 2<=k<=j.
%
%Sintaxis
% [D,err,relerr,n]=diffext(f,x,delta,toler)
%
%Entrada
% - f es la funcion, introducida como una cadena de caracteres 'f'
% - x es el punto en el que se deriva
% - delta es la toleracia para el error
% - toler es la tolerancia para el error relativo
%Resultados:
% - D es la matriz de las aproximaciones de la derivada
% - err es la cota de error
% - relerr es la cota de error relativo
% - n es la coordenada de la "mejor aproximacion"

err=1;
relerr=1;
h=1;
j=1;
D(1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);

while relerr > toler & err > delta &j <12
h=h/2;
D(j+1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
for k=1:j
D(j+1,k+1)=D(j+1,k)+(D(j+1,k)-D(j,k))/((4^k)-1);
end
err=abs(D(j+1,j+1)-D(j,j));
relerr=2*err/(abs(D(j+1,j+1))+abs(D(j,j))+eps);
j=j+1;
end

[n,n]=size(D);