function dfdxdx = num_hesse(f,x,delta=1.0e-3)
## Funktio dfdxdx = num_hesse(f,x,delta) arvioi (symmetrisesti) numeerisesti
## funktion f Hessen matriisin pisteessa x.
##
## Syotteet:
##
## f               Funktio(kahva), jonka Hessen matriisi lasketaan
## x               Piste, jossa Hessen matriisi lasketaan
## delta (1.0e-3)  Pieni luku, jota kaytetaan erotusosamaarassa 
##
## Palautteet:
##
## dfdxdx          Funktion f Hessen matriisi pisteessa x.

if (rows(x) < columns(x)) 
  x=x';
endif

n = length(x);                      # "lyhennysmerkinta"
dfdxdx = zeros(n,n);                # alustus
d = delta*eye(n);                   # differenssimatriisi
for i=1:n
  for j=1:n                         # ALLA PIHVI
    dfdxdx(i,j) = ...
    ( f(x+0.5*d(:,i)+0.5*d(:,j)) - f(x+0.5*d(:,i)-0.5*d(:,j)) ...
     -f(x-0.5*d(:,i)+0.5*d(:,j)) + f(x-0.5*d(:,i)-0.5*d(:,j)) ) / delta^2; 
  endfor
endfor
endfunction