% This function calculates the Hessian matrix of a function ('fun') at point x0 
% This prograpm mimics the procedure "hessp" in the GAUSS library
function y= hessp(fun,x0);
    
% initializations 
[r2, c2] = size(x0);
    k = r2;
    hessian = zeros(k,k);
    grdd = zeros(k,1);
    eps = 6.0554544523933429e-6;

    % Computation of stepsize (dh) for gradient 

    ax0 = abs(x0);
    if x0 ~= 0;
        dax0 = x0./ax0;
    else;
        dax0 = 1;
    end;
   
    
    dh1 = ones(k,2);
    dh1(1:k,1)=ax0;
    dh = eps*max(dh1')'.*dax0;
    xdh = x0+dh;
    dh = xdh-x0;    % This increases precision slightly 
    ee = zeros(k);
    
    i = 1;
    while i <= k;

        ee(i,i) = dh(i);

        i = i+1;
    end;

% Computation of f0=f(x0) 
    f0 = fun(x0);

% Compute forward step 
    i = 1;
    while i <= k;

        grdd(i,1) = fun(x0+ee(:,i));

        i = i+1;
    end;

% Compute "double" forward step 
    i = 1;
    while i <= k;
        j = i;
        while j <= k;

            hessian(i,j) = fun(x0+(ee(:,i)+ee(:,j)));
            if i ~= j;
                hessian(j,i) = hessian(i,j);
            end;

            j = j+1;
        end;
        i = i+1;
    end;

    
grdd1 = zeros(k,k);
f01 = zeros(k,k);
dh1 = zeros(k,k);


i=1;
while i<=k;
        grdd1(:,i)=grdd;
        f01(:,i)=f0;
        dh1(:,i)=dh;
     i=i+1;
end;
    
dh2=dh1.*dh1';
    

y = (((hessian - grdd1) - grdd1') + f0)./dh2 ;