function [hss] = hessjc(f,x0,eps,varargin);

% hessjc. Does hessian. 
% f: function name
% x0: point at which to calculatte
% eps: dx = eps*|x|
%varargin: parameters passed to f
%
if size(x0,2) > 1; 
   disp('Hessjc will bomb -- presumes column vector of input arguments'); 
end; 

f0 = feval(f,x0,varargin{:});
dx = ((x0>0) - (x0<0)).*(eps*x0) + (x0==0).*eps; % if neq0, use dx=e*eps, else use eps
dxv = diag(dx);
hss = zeros(size(x0,1),size(x0,1));
%disp('size(hss)'); disp(size(hss)); 
i = 1;
while i <= size(hss,1);
    %disp('i'); disp(i); 
    hss(i,i) = (feval(f,x0+dxv(:,i),varargin{:})-2*f0+feval(f,x0-dxv(:,i),varargin{:}))/dx(i)^2;
    j = i+1;
    while j <= size(hss,2);
       %disp('j'); disp(j); 
       hss(i,j) =( feval(f,x0+dxv(:,i)+dxv(:,j),varargin{:})-feval(f,x0+dxv(:,i)-dxv(:,j),varargin{:})-...
                 ( feval(f,x0-dxv(:,i)+dxv(:,j),varargin{:})-feval(f,x0-dxv(:,i)-dxv(:,j),varargin{:})))/...
                 (4*dx(i)*dx(j));
              % this uses f + and f -. You can skip some function evaluations by only using f+ and f0. 
        hss(j,i)=hss(i,j);
        j = j+1;
    end;
    i = i+1;
end
return;