function p = Muller(poly,X,Y)

tol=10^(-10);
N=10000;
p1=0;
p2=1/3;
p=2/3;
i=3;
h=inf;
while abs(h)>tol & i<=N
    p0=p1;
    p1=p2;
    p2=p;
    h1=p1-p0;
    h2=p2-p1;
    delta1=(fun(poly,X,Y,p1)-fun(poly,X,Y,p0))/h1;
    delta2=(fun(poly,X,Y,p2)-fun(poly,X,Y,p1))/h2;
    d=(delta2-delta1)/(h2+h1);
    
    b=delta2+h2*d;
    D=(b^2-4*fun(poly,X,Y,p2)*d)^(1/2);
    if abs(b-D)<abs(b+D)
        E=b+d;
    else
        E=b-D;
    end
    h=-2*fun(poly,X,Y,p2)/E;
    p=p2+h;
    i=i+1;
end