function [j1,j2,h1,h2,c1,c2,df1,df2,dj1,dj2,dh1,dh2]=fntest(x,p,perdev1, ...
                                                            perdev2,fname);
% This function checks to see that the jacobian and hessian given as output 
% from fname.m is correctly specified.  It does this by taking numerical  
% derivatives.  I assume that fname.m has inputs (x,p) and outputs 
% (f,jac,hes), the function, jacobian and hessian evaluated at (x,p).
% 'Perdev' should be a vector with length equal to length(x).  The elements 
% of perdev are the percentage deviations away from x used in computing 
% derivatives (e.g. perdev=ones(x)).
%
eval(['[f1,j1,h1,c1,df1,dj1,dh1]=',fname,'(x,p);']);
d=diag(.01*perdev1.*x);
e=diag(.01*perdev2.*p);
arg=find(abs(x)<.0001);
for j=1:length(arg);
  d(arg(j),arg(j))=.001;
end;
arg=find(abs(p)<.0001);
for j=1:length(arg);
  e(arg(j),arg(j))=.001;
end;
nx=length(x);
np=length(p);

for j=1:nx;
  eval(['j2(j,:)=(',fname,'(x+d(:,j),p)-',fname,'(x-d(:,j),p))/(2*d(j,j));']);
  for k=1:j;
    eval(['h2(j,k)=(',fname,'(x+d(:,j)+d(:,k),p)-',fname, ...
      '(x+d(:,j)-d(:,k),p)-',fname,'(x-d(:,j)+d(:,k),p)+', ...
      fname,'(x-d(:,j)-d(:,k),p))/(4*d(j,j)*d(k,k));'])
    h2(k,j)=h2(j,k);
    for l=1:nx;
      eval(['tem=     ',fname,'(x+d(:,j)+d(:,k)+d(:,l),p)-', ...
                        fname,'(x+d(:,j)-d(:,k)+d(:,l),p)-', ...
                        fname,'(x-d(:,j)+d(:,k)+d(:,l),p)+', ...
                        fname,'(x-d(:,j)-d(:,k)+d(:,l),p);']) 
      eval(['tem=tem-(',fname,'(x+d(:,j)+d(:,k)-d(:,l),p)-', ...
                        fname,'(x+d(:,j)-d(:,k)-d(:,l),p)-', ...
                        fname,'(x-d(:,j)+d(:,k)-d(:,l),p)+', ...
                        fname,'(x-d(:,j)-d(:,k)-d(:,l),p));'])
      tem=tem/(8*d(j,j)*d(k,k)*d(l,l)); 
      c2((j-1)*nx+k,l)=tem;
      c2((k-1)*nx+j,l)=tem;
    end;
  end;
end;
if nargout > 6;
dh2=zeros(nx*nx,np);
for j=1:np;      
  eval(['df2(1,j)=(',fname,'(x,p+e(:,j))-',fname,'(x,p-e(:,j)))/(2*e(j,j));']);
  for k=1:nx;
    eval(['dj2(k,j)=(',fname,'(x+d(:,k),p+e(:,j))-',fname, ...
      '(x+d(:,k),p-e(:,j))-',fname,'(x-d(:,k),p+e(:,j))+', ...
      fname,'(x-d(:,k),p-e(:,j)))/(4*d(k,k)*e(j,j));'])
    for l=1:k;
      eval(['tem=     ',fname,'(x+d(:,k)+d(:,l),p+e(:,j))-', ...
                        fname,'(x+d(:,k)-d(:,l),p+e(:,j))-', ...
                        fname,'(x-d(:,k)+d(:,l),p+e(:,j))+', ...
                        fname,'(x-d(:,k)-d(:,l),p+e(:,j));']) 
      eval(['tem=tem-(',fname,'(x+d(:,k)+d(:,l),p-e(:,j))-', ...
                        fname,'(x+d(:,k)-d(:,l),p-e(:,j))-', ...
                        fname,'(x-d(:,k)+d(:,l),p-e(:,j))+', ...
                        fname,'(x-d(:,k)-d(:,l),p-e(:,j)));'])
      tem=tem/(8*d(k,k)*d(l,l)*e(j,j)); 
      dh2((k-1)*nx+l,j)=tem;
      dh2((l-1)*nx+k,j)=tem;
    end;
  end;
end; 
end;
