function x = levmar(F,x0,y,varargin)

    %   x0 is a size N column vector of independent variables to be optimized
    %       to a final value of x
    %   y is a size M column vector of observations to be approached by 
    %       optimizing x
    %   F is a function of N x's which produces M y's, with x and y both
    %       column vectors.  y = feval(F,x,varargin{:})
    %   J is the Jacobian of F w.r.t x, with M rows and N columns.
    %   Therefore J'*J is N by N.  J = jacob(F,x,varargin{:}) ;
    %   ep is a size M column vector which is the error between the observations 
    %       y and feval(F,x,varargin{:}), i.e. ep = y - feval(F,x,varargin{:})
    
    N = length(x0) ; M = length(y) ;
    if M < N, 
        disp('Unable to complete Levenberg-Marquardt algorithm') ;
        disp('Number of observations is less than number of parameters') ;
        return ;
    end
    tau = 1e-3 ; e1 = 1e-15 ; e2 = e1 ; e3 = e1 ;
    k = 0 ; nu = 2 ; x = x0 ; kmax = 200 ;
    J = jacob(F,x,varargin{:}) ;
    A = J'*J ;                          %  A is N by N
    ep = y - feval(F,x,varargin{:}) ;   %  ep is an M by 1 column vector
    g = J'*ep ;                         %  g is an N by 1 column vector
    I = eye(N,N) ;                      %  I is an N by N identity matrix 
    stop = norm(g,inf) <= e1 ; mu = tau*max(trace(A)) ;
    while ~stop & k < kmax,
        k = k + 1 ; rho = -1 ;
        while rho <= 0 & ~stop
            delx = inv(A+mu*I)*g ; %    delx is an N by 1 column vector
            if norm(delx) <= e2*norm(x),
                stop = 1 ;
            else
                xnew = x + delx ;
                rho = (norm(ep)^2 - norm(y-feval(F,xnew,varargin{:}))^2)./...
                        (delx'*(mu*delx+g)) ;
                if rho > 0
                    x = xnew ;
                    A = J'*J ; ep = y - feval(F,x,varargin{:}) ; g = J'*ep ;
                    stop = (norm(g,inf) <= e1) | (norm(ep)^2 <= e3) ;
                    mu = mu*max(1/3,1-(2*rho-1)^3) ; nu = 2 ;
                else
                    mu = mu*nu ; nu = 2*nu ;
                end
            end
        end
     end
 
 function J = jacob(F,x,varargin)
    y = feval(F,x,varargin{:}) ; M = length(y) ; N = length(x) ; J = zeros(M,N) ;
    for n = 1:N,
        x2 = x ; dx = max(abs(x(n))/1000,1e-6) ; 
        x2(n) = x(n) + dx ; y2 = feval(F,x2,varargin{:}) ;
        J(:,n) = (y2 - y)/dx ;
    end