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