function dfdx = num_jacobi(f,x,delta=1.0e-3)
## Funktio dfdx = num_jacobi(f,x,delta) arvioi (symmetrisesti) numeerisesti
## funktion f Jacobin matriisin (osittaisderivaattamatriisin) 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 Jacobin matriisi pisteessa x.

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