function [z,n,delta] = muller(p,q,r,lam,nmax,tol)
n = 0;z = lam;p2 = 2*p;delta = 1;
nmax
tol
while(n <= nmax & abs(delta) > tol*z)
    n = n+1;
    z2  = z*z;
    
    fp  = (4*z2 + p2)*z + q;
    f   = z2*(z2 + p) + q*z + r;
    if (n <=10)   %take a simple Newton step to improve guess
        delta = -f/fp; 
    else         %proceed with quadratic corrections (not worth it!)
        fpp = 12*z2 + p2;
        delta = (-fp + (fp*fp - 2*fpp*f)^(1/2))/fpp;
    end
    z = z + delta
end
% Newton takes two extra steps to converge---at less cost