function y=muller(f,x0,x1,x2,delta,maxsteps)

% y=MULLER(f,x0,x1,x2,delta,maxsteps) dela Mullerjevo metodo za funkcijo f 
% z zacetnimi priblizki x0, x1 in x2. Iteracija se konca, ko se zadnja dva 
% priblizka razlikujeta za manj kot delta ali ko prekoracimo maksimalno 
% stevilo korakov maxsteps. 
%
% Primer uporabe:
% f = @(x) x^3-5*x+1
% muller(f,2,2.05,2.1,1e-12,50)

% Bor Plestenjak 2020

% privzete vrednosti argumentov, če jih ne podamo
if nargin<5, delta=1e-12; end
if nargin<6, maxsteps=50; end

xn=x2;  x2=x1;  x1=x0;  korak=0;                    						 

while (abs(xn-x2)>delta) && (korak<maxsteps)
    x0 = x1;
    x1 = x2;
    x2 = xn;
    ulomek = (x0-x2)*(x1-x2)*(x0-x1);
    a = ( (x1-x2)*(f(x0)-f(x2))-(x0-x2)*(f(x1)-f(x2)) )/ulomek;
    b = ( (x0-x2)^2*(f(x1)-f(x2))-(x1-x2)^2*(f(x0)-f(x2)) )/ulomek;
    c = f(x2);
    D = sqrt(b^2-4*a*c);
    if b<0 
        xn = x2 - 2*c/(b-sqrt(b^2-4*a*c));
    else
        xn = x2 - 2*c/(b+sqrt(b^2-4*a*c));
    end        
    korak = korak+1;
    if (imag(xn)~=0)
       disp(sprintf('%3d: %15.15e %+15.15ei  %0.1e',[korak,real(xn),imag(xn),abs(xn-x2)]));  
    else
       disp(sprintf('%3d: %15.15e  %0.1e',[korak,xn,abs(xn-x2)]));  
    end
    pause(0.2)
end   

y=xn;