function p=newton(f,df,p0,tol)

while 1
  p=p0-f(p0)/df(p0);
  if abs(p-p0)<tol, break; end
  p0=p;
end
