function dfdx = num_grad(f,x,delta=1.0e-3)
## Funktio dfdx = num_grad(f,x,delta) arvioi (symmetrisesti) numeerisesti
## funktion f gradientin pisteessa x.
##
## Syotteet:
##
## f               Funktio(kahva), jonka gradientti lasketaan
## x               Piste, jossa gradientti lasketaan
## delta (1.0e-3)  Pieni luku, jota kaytetaan erotusosamaarassa 
##
## Palautteet:
##
## dfdx            Funktion f gradientti pisteessa x.

if (rows(x) < columns(x)) 
  x=x';                                          # varmistetaan pystyvektorius
endif
n = length(x);                                   # "lyhennysmerkinta"
dfdx = zeros(n,1);                               # alustus
d = delta*eye(n);                                # differenssimatriisidiagonaali
for i=1:n
  dfdx(i) = ...
  ( f(x+0.5*d(:,i)) - f(x-0.5*d(:,i)) ) / delta; # PIHVIN PAIKKA
endfor
endfunction