function [p,yp,err,Q] = steff(f,df,p0,delta,epsilon,max1)
% [p,yp,err] = steff(f,df,p0,delta,epsilon,max1)
% [p,yp,err,Q] = steff(f,df,p0,delta,epsilon,max1)
% Steffensen's method is used to locate a root.
% f is the function, input.
% d is the derivative of f, input.
% p0 is the starting point, input.
% delta is the tolerance for p, input.
% epsilon is the tolerance for yp, input.
% max1 is the maximum number of iterations, input.
% p is the root, output.
% yp is the function value, output.
% err is the error estimate for p, output.
% Q is the is the vector of iterations, output.
y0 = feval(f,p0);
p = p0; yp = y0;
Q(1) = 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);
p = p1; yp = y1;
Q = [Q,p1];
err = abs(dp);
relerr = err/(abs(p1)+eps);
if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
df1 = feval(df,p1);
if df1 == 0,
dp = 0;
else
dp = y1/df1;
end
p2 = p1 - dp;
y2 = feval(f,p2);
p = p2; yp = y2;
Q = [Q,p2];
err = abs(dp);
relerr = err/(abs(p1)+eps);
if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
d1 = (p1 - p0)^2;
d2 = p2 - 2*p1 + p0;
if d2 ==0,
dp = p2 - p1;
p3 = p2;
break;
else
dp = d1/d2;
p3 = p0 - dp;
end
p0 = p3;
y0 = feval(f,p0);
p = p0; yp = y0;
Q = [Q,p0];
err = abs(dp);
relerr = err/(abs(p3)+eps);
if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
end