function [b,q]=numder2(f,x,par,perdev)
%NUMDER calculates the first and second two-sided numerical derivatives of 
%  a function f(x,par) with respect x.
%
%  INPUTS:
%    (i) f= a string giving filename for function.  The first output argument
%           of `f'.m should be a scalar giving f(x)
%   (ii) x= argument of f(.) and point at which derivatives evaluated.
%  (iii) par= parameter vector needed by f(.).  If none required, set up `f'.m
%             with dummy argument for parameters and set par=0.
%   (iv) perdev= a vector (of same length as x) giving percent deviations
%                to take when computing derivatives  
%
%  OUTPUTS:
%    (i) b= df(x)/dx, the numerical first derivative of f(x) at x.
%   (ii) q=d/dx'(df(x)/dx), the numerical second derivative of f(x) at x.
%
d=diag(.01*perdev.*x);
arg=find(abs(x)<1e-9);
for j=1:length(arg);
  d(arg(j),arg(j))=1e-4;
end;
n=length(x);
for j=1:n;
  eval(['fxp(j)=',f,'(x+d(:,j),par);']);
  eval(['fxm(j)=',f,'(x-d(:,j),par);']);
  b(j,1)=(fxp(j)-fxm(j))/(2*d(j,j));
end;
if nargout > 1;
  for j=1:n;
    eval(['q(j,j)=(fxp(j)-2*',f,'(x,par)+fxm(j))/d(j,j)^2;']);
    for m=1:j-1;
        eval(['q(j,m)=(',f,'(x+d(:,j)+d(:,m),par)-',f,'(x+d(:,j)-d(:,m),', ..
            'par)-',f,'(x-d(:,j)+d(:,m),par)+',f,'(x-d(:,j)-d(:,m),par))/',  ..
            '(4*d(j,j)*d(m,m));']);
    end;
  end;
  q=q+tril(q,-1)';
end;
