function [yj1,yj2]=fncheck(x,param,perdev,fname);
% This function checks to see that the jacobian given as the second output
% from fname.m is correctly specified. It does this by taking numerical
% derivatives. I assume that fname.m has inputs (x,param) and outputs
% (f,fjac,fhes), where the function and jacobian are evaluated at (x,param).
% 'Perdev' should be a matrix equal in size x. The elements
% of perdev are the percentage deviations away from x used in computing
% derivatives (e.g. perdev=ones(x)).
%
eval(['[y1,yj1]=',fname,'(x,param);']);
[m,n]=size(x);
z=zeros(m,n);
for j=1:n;
for i=1:m;
dev=z;
dev(i,j)=max(.01*perdev(i,j)*x(i,j),1e-8);
eval(['term1=',fname,'(x+dev,param);'])
eval(['term2=',fname,'(x-dev,param);'])
yj2(:,(j-1)*m+i)=(term1(:)-term2(:))/(2*dev(i,j));
end;
end;