Matlab - metodo extrapolacion de richardson para derivacion numerica

Hola! me preguntaria si alguien podria hacerme el favor de ayudarme a hacer un programa en matlab para aproximar numericamente la matriz jacobiana (formada por las derivadas parciales de primer orden de una funcion ) mediante el metodo de extrapolacion de richardson ya q eh intentado hacerlo pero es bastante difícil por que es la primera vez que estoy viendo programacion en analisis numerico.
Gracias por su colaboracion.

aqui hay una ayuda pa q se guien de un programa donde se aproxima numericamente la derivada de cualquier funcion en una sola variable mediante la extrapolacion de richardson.

function [D,err,relerr]=richardson(f,x0,h,tol)
%% entradas
% f es la funcion
% x0 es el punto en el q se deriva
% h es la tolerancia para el error
% tol es la tolerancia para el error relativo

%% salidas
% D es la matriz de las aproximaciones a las derivada
% err es la cota del error
% relerr es la cota del error relativo
format long
err=1;
relerr=1;
j=1;
x=x0+h;
f1=eval(f);
x=x0-h;
f2=eval(f);
D(1,1)=(f1-f2)/(2*h);
while (err>h & relerr>tol)
h=h/2;
x=x0+h;
f3=eval(f);
x=x0-h;
f4=eval(f);
D(j+1,1)=(f3-f4)/(2*h);
for k=1:j
D(j+1,k+1)=((4^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);