function y=muller(p0,p1,p2)
%
% x0 = muller(p0,p1,p2) USES MULLER'S METHOD TO
% APPROXIMATE THE ROOT OF THE FUNCTION f(x) GIVEN
% THE NUMBERS p0, p1 AND p2.  THE ANSWER IS ACCURATE
% TO WITHIN tol
%

% STOPPING CRITERIA
tol = 1e-8; N0 = 50;

% INITIAL POINTS
y(1) = p0;   y(2) = p1;   y(3) = p2;

% THE FUNCTION
Fmul =@(x) x.^4 - 3*x.^3+x.^2 + x + 1;

% START THE ITERATIONS USING THE BOOK'S ALGORITHM
i = 3; 

while i <= N0,
   C = Fmul(p2);
   h1 = p1-p0;
   h2 = p2-p1;
   d1 = (Fmul(p1)-Fmul(p0))/h1;
   d2 = (C-Fmul(p1))/h2;
   A = (d2-d1)/(h2+h1);
   B = d2+h2*A;
   Disc = sqrt(B^2-4*A*C);
   if abs(B-Disc)<abs(B+Disc), E = B+Disc; else E = B-Disc; end
   h = -2*C/E;

   p = p2 + h;
   y(i+1) = p;
   if abs(h) < tol,  break;   end
   i = i + 1;
   p0 = p1;
   p1 = p2; 
   p2 = p; 
   h1 = h2;
   h2 = h;
end
y=y(:);

% DISPLAY THE SOLUTION ON THE SCREEN
format long
if i<=N0,
   disp(' ')
   disp(['The root is ',num2str(y(end),9),' with a possible error of ',num2str(tol)])
   disp(' ')
   disp(' n                   p_n                     |F(p_n)|         ')       
   disp('--------------------------------------------------------------')
   
   for kk=1:length(y)
      disp(sprintf('%2.0f   %15.12f + %15.12fi   %15.12e   ',kk,[real(y(kk)) imag(y(kk))],abs(Fmul(y(kk))) )) 
   end
   disp(' ')
else
   disp(' ')
   disp(['Muller''s method failed after ',num2str(N0),' iterations']),
   disp(' ')
end
