function [hes_x] = hessian(fun_name, nor_x, error)

header;

% [hes_x] = hessian(fun_name, nor_x, error)
% This function finds the numerical hessian of a function at a point
% 'hes_x' is the hessian returned by the function
% 'fun_name' is the function for which we want to find the hessian
% 'nor_x' is the point where we want to find the hessian
% 'error' is the demanded accuracy
%
% File hessian.m
% ESCI 416, Term Project, Fall 1995, by Kittipan Techakittiroj and Jonas Karlsson

N = length(nor_x);
variation = zeros(size(nor_x));
hes_x = zeros(length(nor_x), length(nor_x));

for i = 1:N,
  variation(i) = error;
  before = gradient(fun_name, nor_x - variation, error);
                                                    gra_eval = gra_eval + 1;
                                                    add_vector = add_vector + 1;
  after = gradient(fun_name, nor_x + variation, error);
                                                    gra_eval = gra_eval + 1;
                                                    add_vector = add_vector + 1;
  hes_x(i,:) = (after - before) / (2 * error);
                                                    mul_scalar = mul_scalar + 1;
                                                    div_scalar = div_scalar + 1;
                                                    add_vector = add_vector + 1;
                                                    c_m_vector = c_m_vector + 1;
  variation(i) = 0;
end;