function [p0,y0,err,P] = newton(f,df,p0,delta,epsilon,max1)
% [p0,y0,err] = newton(f,df,p0,delta,epsilon,max1)
% [p0,y0,err,P] = newton(f,df,p0,delta,epsilon,max1)
% Newton's method is used to locate a root.
% f is the function, input.
% df is the derivative of f, input.
% p0 is the starting point, input.
% delta is the tolerance for p0, input.
% epsilon is the tolerance for y0, input.
% max1 is the maximum number of iterations, input.
% p0 is the root, output.
% y0 is the function value, output.
% err is the error estimate for p0, output.
% P is the is the vector of iterations, output.
P(1) = p0;
y0 = feval(f,p0);
for k=1:max1,
df0 = feval(df,p0);
if df0 == 0,
dp = 0;
else
dp = y0/df0;
end
p1 = p0 - dp;
y1 = feval(f,p1);
err = abs(dp);
relerr = err/(abs(p1)+eps);
p0 = p1;
y0 = y1;
P = [P;p1];
if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
end