function [J, R, V, template] = jacobian(A);% Use finite-differences to compute Jacobian% Copyright (c) by Raymond A. Adomaitis, 1998-2005Akey = keys(get(A,'var'));Aval = values(get(A,'var'));[V, template] = columnnae(A,'var');Aresjacflag = get(A,'resjacflag');% Finite difference scaling is based on Solving Nonlinear Equations% with Newton's Method, C. T. Kelly, SIAM, 2003.Epsil      = sqrt(eps)*max([abs(V),ones(size(V))],[],2).*sgn(V);VPlusEpsil = V + Epsil;if ~Aresjacflag     % default - residual is not set up to                    % compute the Jacobian by finite differences% check if residual in object A must be updated   if ~get(A,'residflag'), A = residual(A); end   R = columnnae(A,'resid');   n = 0;   B = A;   for i = 1:length(Akey)      for j = 1:prod(template{i})         n = n+1;         B = set(B,'var', ...                values(setval(get(A,'var'),VPlusEpsil(n),Akey{i},j)));         B = residual(B);                     RPlusEpsil = columnnae(B,'resid');         J(:,n) = ( RPlusEpsil - R )/Epsil(n);               end      endelse % residual is set up to compute Jacobian by finite differences   [A, J, R] = residual(A,Akey,Aval,V,template, Epsil,VPlusEpsil);   end% RAA 9/15/2004 changed from: %%    Epsil = 1e-4*mean(abs(V));%    if Epsil < sqrt(eps), Epsil = sqrt(eps); end%    VDiff = Epsil*V;%    I = find( abs(VDiff) < Epsil );%    VPlusEpsil(I) = V(I) + Epsil;%    VDiff(I) = Epsil;%    J(:,n) = ( RPlusEpsil - R )/VDiff(n);% -------------------- sgn.m -------------------- %function B = sgn(A)% return the sign of A if nonzero; return 1 if zeroI = find( A==0 );if ~isempty(I)   B(I,1) = 1;   endJ = find( A~=0 );if ~isempty(J)   B(J,1) = A(J,1)./abs(A(J,1));   end