function [p2,y2,err,P] = muller(f,p0,p1,p2,delta,epsilon,max1)
% [p2,y2,err] = muller(f,p0,p1,p2,delta,epsilon,max1)
% [p2,y2,err,P] = muller(f,p0,p1,p2,delta,epsilon,max1)
% Muller's method is used to locate a root.
% f is the function, input.
% p0 is a starting point, input.
% p1 is a starting point, input.
% p2 is a starting point, input.
% delta is the tolerance for p2, input.
% epsilon is the tolerance for y2, input.
% max1 is the maximum number of iterations, input.
% p2 is the root, output.
% y2 is the function value, output.
% err is the error estimate for p2, output.
% P is the is the vector of iterations, output.
P(1) = p0;
P(2) = p1;
P(3) = p2;
y0 = feval(f,p0);
y1 = feval(f,p1);
y2 = feval(f,p2);
for k=1:max1,
h0 = p0 - p2;
h1 = p1 - p2;
c = y2;
e0 = y0 - c;
e1 = y1 - c;
det1 = h0*h1*(h0-h1);
a = (e0*h1 - h0*e1)/det1;
b = (h0^2*e1 - h1^2*e0)/det1;
if b^2 > 4*a*c,
disc = sqrt(b^2 - 4*a*c);
else
disc = 0;
end
if b < 0, disc = - disc; end
z = - 2*c/(b + disc);
p3 = p2 + z;
if abs(p3-p1) < abs(p3-p0),
u = p1;
p1 = p0;
p0 = u;
v = y1;
y1 = y0;
y0 = v;
end
if abs(p3-p2) < abs(p3-p1),
u = p2;
p2 = p1;
p1 = u;
v = y2;
y2 = y1;
y1 = v;
end
p2 = p3;
y2 = feval(f,p2);
P = [P,p2];
err = abs(z);
relerr = err/(abs(p3)+eps);
if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
end