function [g,H] = Hessian(f,x,h,varargin)
% PURPOSE: Computes finite difference Hessian
% -------------------------------------------------------
% Usage:  H = hessian(func,x,varargin)
% Where: func = function name, fval = func(x,varargin)
%           x = vector of parameters (n x 1)
%    varargin = optional arguments passed to the function
% -------------------------------------------------------
% RETURNS:
%           H = finite differnce hessian
% -------------------------------------------------------

eps = 1e-5;
n = size(x,1);
%h=1e-5;
xh = x+h;
h = xh-x;    
ee = sparse(1:n,1:n,h,n,n);
% Central Difference Approximations
f_counter = 0;
% First-order derivatives: 2n function calls needed

%for i=1:n
%    f_counter = f_counter + 1;
%    f_evaluate(f_counter,:) =  x+ee(:,i);
%    %A = feval(f,x+ee(:,i),varargin{:});
%    
%    f_counter = f_counter + 1;
%    f_evaluate(f_counter,:) =  x-ee(:,i);
%    %B = feval(f,x-ee(:,i),varargin{:});
%end

for i=1:n
    A = feval(f,x+ee(:,i),varargin{:});
    A = sum(A.^2)/size(A,1)
    B = feval(f,x-ee(:,i),varargin{:});
    B = sum(B.^2)/size(B,1)
    g(i) = (A-B)/(2*ee(i,i))
end
%Second-order derivatives based on function calls only (Abramowitz and Stegun 1972, p.884): for dense Hessian, 2n+4n^2/2 function calls needed.

%for i=1:n
%    f_counter = f_counter + 1;
%    f_evaluate(f_counter,:) =x+2*ee(:,i);
%    %C = feval(f,x+2*ee(:,i),varargin{:});
%    
%    f_counter = f_counter + 1;
%    f_evaluate(f_counter,:) =x;
%    %E = feval(f,x,varargin{:});
%    
%    f_counter = f_counter + 1;
%    f_evaluate(f_counter,:) =  x-2*ee(:,i);    
%    %F = feval(f,x-2*ee(:,i),varargin{:});
%    for j=i:n
%        f_counter = f_counter + 1;
%        f_evaluate(f_counter,:) =x+ee(:,i)+ee(:,j);
%        %G = feval(f,x+ee(:,i)+ee(:,j),varargin{:});
%        
%        f_counter = f_counter + 1;
%        f_evaluate(f_counter,:) =x+ee(:,i)-ee(:,j);            
%        %I = feval(f,x+ee(:,i)-ee(:,j),varargin{:});
%        
%        f_counter = f_counter + 1;
%        f_evaluate(f_counter,:) =x-ee(:,i)+ee(:,j);            
%        %J = feval(f,x-ee(:,i)+ee(:,j),varargin{:});
%        
%        f_counter = f_counter + 1;
%        f_evaluate(f_counter,:) =x-ee(:,i)-ee(:,j);            
%        %K = feval(f,x-ee(:,i)-ee(:,j),varargin{:});
%    end
%end

for i=1:n
    C = feval(f,x+2*ee(:,i),varargin{:});
    C = sum(C.^2)/size(C,1);
    E = feval(f,x,varargin{:});
    E = sum(E.^2)/size(E,1);
    F = feval(f,x-2*ee(:,i),varargin{:});
    F = sum(F.^2)/size(F,1); 
    H(i,i) = (- C + 16*A - 30*E + 16*B - F)/(12*(ee(i,i)^2));
    for j=i:n
        G = feval(f,x+ee(:,i)+ee(:,j),varargin{:});
        G = sum(G.^2)/size(G,1)   
        I = feval(f,x+ee(:,i)-ee(:,j),varargin{:});
        I = sum(I.^2)/size(I,1)   
        J = feval(f,x-ee(:,i)+ee(:,j),varargin{:});
        J = sum(J.^2)/size(J,1)           
        K = feval(f,x-ee(:,i)-ee(:,j),varargin{:});
        K = sum(K.^2)/size(K,1)    
        H(i,j) = (G-I-J+K)/(4*ee(i,i)*ee(j,j));
        H(j,i) = H(i,j);
    end
end
%disp(H)
%Second-order derivatives based on gradient: 2n gradient calls needed.
%disp(f_counter)
%disp(f_evaluate)