function [x,polku] = ...
newton_stat(f,x0,polkuko=0,maxit=10,hyppy=1.0,delta=1.0e-3)
## Funktio [x,polku] = newton_sta(f,x0,polkuko=0,maxit=10,hyppy=1.0,delta=1.0e-3) 
## etsii reaaliarvoisen funktion gradientin nollakohdan (eli stationaarisen 
## pisteen) Newtonin menetelmalla.
##
## Syotteet:
##    
## f               Funktio(kahva) f, jonka maksimia etsitaan. 
## x0              Alkuarvaus maksipisteelle.
## polkuko    (0)  Palautetaan koko iteraatiopolku, jos =1 (oletusarvo =0). 
## maxit     (10)  Iteraatioiden lukumaara.
## hyppy    (1.0)  Iteraatohypyn (suhteellinen) koko.
## delta (1.0e-3)  Numeerisessa derivoinnissa kaytettava differenssi.
##
## Palautteet:
##
## x               Arvio funktion f gradientin nollakohdasta.
## polku           Iteraatiopolku, jos polkuko=1, muulloin x0.

if (rows(x0) < columns(x0)) 
  x0=x0';                                # varmistetaan pystyvektorius
endif

polku=x0;
if (polkuko==1)
  polku = rand(length(x0),maxit);        # alustetaan polku (tarvittaessa)
endif

x_vanha = x0;
for i=1:maxit
  if (polkuko==1)
    polku(:,i)=x_vanha;                  # kerataan polku talteen (tarvittaessa)
  endif
  H = num_hesse(f,x_vanha,delta);        # Hesse numeerisesti
  grad = num_grad(f,x_vanha,delta);      # gradientti numeerisesti

  x_uusi = x_vanha - hyppy*inv(H)*grad;  # PIHVIN PAIKKA
  x_vanha = x_uusi;
endfor
x = x_uusi;
endfunction