function [Jac] = zgn_numdif( fcn, x, xthrsh, uds, F)
%
%  M-function zgn_iniscalx performs numerical differentiation
%  to get the Jacobian J(x) for NLSQ_RES and NLSQ_ERR
%
%
  n = length( x); 
  xh = x;
  for icol = 1:n;
    delx = max( abs(x(icol)), xthrsh( icol)) * 1.e-7;
    xh(icol) = x(icol) + delx;
    Fh = feval( fcn, xh, '', uds);
    Jac(:,icol) = (Fh - F) / delx;
    xh(icol) = x(icol);
  end
