function p=muller(f,p0,p1,p2,tol,max_iterations)

% Mueller's method: See section 2.6.

% Variable names:
% p: approximate value of root
% h: difference from previous approximation
% d: slope of P(x) between previous and current p-values
% a, b, c: quadratic coefficients (in terms of p-h)
h1=p1-p0;
h2=p2-p1;
d1=(f(p1)-f(p0))/h1;
d2=(f(p2)-f(p1))/h2;
for n=1:max_iterations
    a=(d2-d1)/(h2+h1);
    b=d2+h2*a;
    c=f(p2);
    % h = quadratic formula
    D=sqrt(b^2-4*a*c);
    if abs(b-D)<abs(b+D)
        E=b+D;
    else
        E=b-D;
    end
    h=-2*c/E;
    % done with quadratic formula
    p=p2+h;
    if abs(p-p2)<tol
        return
    end
    p0=p1;
    p1=p2;
    p2=p;
    h1=p1-p0;
    h2=p2-p1;
    d1=(f(p1)-f(p0))/h1;
    d2=(f(p2)-f(p1))/h2;
end