eval(['[L,dL]=',func,'(x);'])
d=diag(.0001*abs(x));
n=length(x);
for j=1:n;
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  J(j,1)=(fxp-fxm)/(2*d(j,j));
end;

y=find(abs(J-dL)./abs(dL) > 1e-3);
m=length(y);

for k=1:m;
  j=y(k);
  d=diag(.1*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,1)=(fxp-fxm)/(2*d(j,j));
  d=diag(.01*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,2)=(fxp-fxm)/(2*d(j,j));
  d=diag(.001*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,3)=(fxp-fxm)/(2*d(j,j));
  d=diag(.0001*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,4)=(fxp-fxm)/(2*d(j,j));
  d=diag(.00001*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,5)=(fxp-fxm)/(2*d(j,j));
  d=diag(.000001*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,6)=(fxp-fxm)/(2*d(j,j));
  d=diag(.0000001*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,7)=(fxp-fxm)/(2*d(j,j));
  d=diag(.00000001*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,8)=(fxp-fxm)/(2*d(j,j));
  d=diag(.000000001*abs(x)); 
  eval(['fxp=',func,'(x+d(:,j),par);']);
  eval(['fxm=',func,'(x-d(:,j),par);']);
  JJ(k,9)=(fxp-fxm)/(2*d(j,j));
end;
