function dm = demo_muller()

    eps = 1.0e-6;

    x0 = 1.0;
    x1 = 2.0;
    x2 = 1.1;

    y0 = poly(x0);
    y1 = poly(x1);
    y2 = poly(x2);

    n = 1;
    while ( abs(y2) > eps ) 
    
        fprintf('Krok %3d:  koren x= %8.4e,  f(x)= %8.4e\n', n, x2 ,y2);

        a0 = y0 / ( (x0-x1)*(x0-x2) );
        a1 = y1 / ( (x1-x0)*(x1-x2) );
        a2 = y2 / ( (x2-x0)*(x2-x1) );

        A = a0+a1+a2;
        C = x1*x2*a0 + x0*x2*a1 + x0*x1*a2;
        B = -(x1+x2)*a0 -(x0+x2)*a1 -(x0+x1)*a2;

        xn   = ( -B + sqrt(B^2-4*A*C) ) / (2*A);
        xn0  = ( -B - sqrt(B^2-4*A*C) ) / (2*A);

        if ( abs(xn0-x2) < abs(xn-x2) ) 
            xn = xn0;
        end;
            
        x0 = x1;
        x1 = x2;
        x2 = xn;
        y0 = y1;
        y1 = y2;
        y2 = poly(xn);

        n = n+1;
        pause;
    end;

    return;
    

%%%%%%%%%%%%

function p = poly(x);

    p = 4*x^3 - 2*x^2 - 4*x - 3;

    return;

%%%%%%%%%%%%
