function A = newton(A,paramlist)% newton.m A Newton-Raphson procedure for computing solutions%            to naemodel-derived objects. The method uses a%            finite-difference approximation for the Jacobian%            array (see the jacobian.m method) and updates both%            the 'var' and 'resid' fields of the neqmodel-derived%            object. Called as%%            A = newton(A,paramlist)%% INPUT/OUTPUT PARAMETERS%            A : a naemodel-derived object%    paramlist : a cell-array of parameter names from A.param;%                if specified, sensitivity.m is used instead of%                jacobian.m% Copyright (c) by Raymond A. Adomaitis, 1998-2004unpack(A,'solverset');if isempty(vmin) & isempty(vmax)    checklim = 0;else    checklim = 1;    vmin = columnnae(A,'solver_vmin');    vmax = columnnae(A,'solver_vmax');endif lsmax < 2    warning('Minimum value for lsmax = 2; reset to 2 in this object')    A = set(A,'solverset',2,'lsmax');    lsmax = get(A,'solverset','lsmax');endfor inewt = 1:nmax    % note that the jacobian method updates the residual in object A    if inewt == 1 | ~mod(inewt,jup) % Newton (if true), modified Newt (false)        if nargin == 2            [J, R, param, template] = sensitivity(A,paramlist);            sjflag = 1;        else            [J, R, var, template] = jacobian(A);            sjflag = 0;        end        % 'normalize' each column of J        for i = 1:size(J,2)            colnorm(i) = mean(abs(J(:,i)));            if colnorm(i) < eps                error(['Column ',int2str(i),' of ',int2str(size(J,2)), ...                    ' in the Jacobian has vanished'])            else                J(:,i) = J(:,i)/colnorm(i);            end        end    else % do not update Jacobian at this step        A   = residual(A);        R   = columnnae(A,'resid');        var = columnnae(A,'var');    end        if size(J,1) > size(J,2)        update = -J\R;    else        if rcond(J) > sqrt(eps)            update = -J\R;        else            update = linearsystemsolver(-J,R);        end    end        if ~(abs(update) >= 0), error('Exiting the newton method'), end    update = real( update./colnorm' );    % Perform a line search in the direction of the update vector    % to determine approximate minimum residual value    Rnorm(1) = norm(R);    uplength(1) = 0;    uplength(2) = 1.0;    for i = 2:lsmax        if sjflag            paramup = col2cell(param + uplength(i)*update,template);            B = set(A,'param',paramup);        else            if checklim                % Check is current uplength will generate an update that exceeds                % any specified limits on variable values                % find indices of updated variables that exceed lower bound                exlb = find(var + uplength(i)*update <= vmin);                if ~isempty(exlb)                    frac = (vmin(exlb)-var(exlb))./(uplength(i)*update(exlb));                    [minfrac ifrac] = min(abs(frac));                    uplengthlimit = sign(frac(ifrac))*minfrac*uplength(i);                    uplength(i) = (uplengthlimit + uplength(i-1))/2;                end                % now check upper bound                exub = find(var + uplength(i)*update >= vmax);                if ~isempty(exub)                    frac = (vmax(exub)-var(exub))./(uplength(i)*update(exub));                    [minfrac ifrac] = min(abs(frac));                    uplengthlimit = sign(frac(ifrac))*minfrac*uplength(i);                    uplength(i) = (uplengthlimit + uplength(i-1))/2;                end            end % ends "if checklim"            varup = col2cell(var + uplength(i)*update,template);            B = set(A,'var',varup);        end % ends "if sjflag"        B = residual(B);        R = columnnae(B,'resid');        Rnorm(i) = norm(R);        if Rnorm(i) <= Rnorm(i-1)            uplength(i+1) = uplength(i) + 2/3*(uplength(i)-uplength(i-1));        elseif Rnorm(i) > Rnorm(i-1)            uplength(i+1) = uplength(i) - 2/3*(uplength(i)-uplength(i-1));        else            error('Rnorm is NaN - quitting Newton')        end    end % ends i loop    A = B;    if get(A,'solverset','plot')        plot(A), drawnow    end    disp([ 'Update/residual norm = ', ...        num2str(norm(uplength(lsmax)*update)),' / ', ...        num2str(Rnorm(lsmax)) ])    if inewt == 1, normup1 = norm(uplength(lsmax)*update); end    if norm(uplength(lsmax)*update)/normup1 < 1e-12 | ...            norm(uplength(lsmax)*update) < 1e-12        break    endend