function [p0, err, i, y] = newton(f,df,p0,toly,tolx,max)
for i=1:max
   p1=p0-feval(f,p0)/feval(df,p0);
   err=abs(p1-p0);
   relerr=2*err/(abs(p1)+eps);
   p0=p1;
   y=feval(f,p0);
   if (err<tolx) | (relerr<tolx) | abs(y)<toly, break, end
end
