function  [J,H]=numder(f,x,par,perdev,side)
%NUMDER   Calculates numerical first and second derivatives.
%         [J,H]=numder(f,x,par,perdev,side) calculates the Jacobian and
%         hessian for the scalar-valued function f:
%
%         J=[ df/dx1 ]   H=[ d/dx1(df/dx1) d/dx1(df/dx2) .. d/dx1(df/dxn) ]
%           [ df/dx2 ]     [ d/dx2(df/dx1) d/dx2(df/dx2) .. d/dx2(df/dxn) ]
%           [   :    ]     [      :                              :        ]
%           [   :    ]     [      :                              :        ]
%           [ df/dxn ]     [ d/dxn(df/dx1) d/dxn(df/dx2) .. d/dxn(df/dxn) ]
% 
%         where x=[x1, x2, x3, ... xn]' and f has input arguments x and
%         par.   The argument perdev is a n x 1 vector where perdev(i) is 
%         the percent deviation for xi used when computing difference 
%         intervals.  The argument s can be a scalar equal to 1 or 2 or a 
%         n x 1 vector with s(i) = 1 or 2, i=1,...n.  If s(i)=1, one-sided
%         derivatives are computed and if s(i)=2, two-sided derivatives
%         are computed for xi.
%
%         [J,H]=numder(f,x,par) uses perdev(i) = max(.0001*abs(x(i)), 1e-8)
%         for all i and two-sided derivatives.
%

%         Ellen R. McGrattan  4-5-87
%         Revised, 3-2-90, ERM

x  = x(:);
n  = length(x);
if nargin <4;
  side  = 2*ones(n,1);
  d     = diag(max (.0001*abs(x),1e-8*ones(n,1)));
end; 
if side==1 & n>1; side =   ones(n,1); end;
if side==2 & n>1; side = 2*ones(n,1); end;
d  = diag(max(.01*perdev.*abs(x),1e-8*ones(n,1)));
sd = diag(side-1)*d;
for j=1:n;
  eval(['fxp(j)=',f,'(x+d(:,j),par);']);
  eval(['fxm(j)=',f,'(x-sd(:,j),par);']);
  J(j,1) = (fxp(j)-fxm(j))/(side(j)*d(j,j));
end;
if nargout > 1;
  H=zeros(n);
  for i=1:n;
    for j=1:i;
      eval(['H(i,j)=(',f,'(x+d(:,i)+d(:,j),par)-',f,'(x+d(:,i)-sd(:,j),', ...
           'par)-',f,'(x-sd(:,i)+d(:,j),par)+',f,'(x-sd(:,i)-sd(:,j),par))', ...
           '/(side(i)*side(j)*d(i,i)*d(j,j));']);
    end;
  end;
  H=H+tril(H,-1)';
end;
