function [J]=Jac(x,At,b) % This function modified from a code written by C. Fuhrer FMN081-2005
% computes the Jacobian of a function 
n=length(x); 
fx=F(x,At,b); 
eps=1.e-8; % could be made more accurate
xperturb=x; 
for i=1:n 
    xperturb(i)=xperturb(i)+eps; 
    J(:,i)=(F(xperturb,At,b)-fx)/eps; 
    xperturb(i)=x(i); 
end;