% Root finding using MULLER METHOD. Modified Sep 28, 2004

function Muller (Poly, x1,x2,x3,A, B,Tol,itmax)

fprintf(' MULLER Method.\n');
format short e

% Input data.

disp('   Input data    ');
disp('        x1        x2             x3           Tol       itmax');
disp( [x1,x2,x3,Tol, itmax] );

% Graph of the Polynomial for which the roots are being found.

step=(B-A)/100;
x=A:step:B;
y=feval(Poly,x);
plot(x,y)
xlabel('x'); ylabel('f(x)'); 
grid on

% **************************MAIN   LOOP ****************************
 it=0;
   while it < itmax 
       it = it + 1;
       h1=x1-x2;
       h2=x2-x3;
       h3=x1-x3;
       Px1=feval(Poly,x1);
       Px2=feval(Poly,x2);
       Px3=feval(Poly,x3);
       delta2=Px2-Px3;
       delta3=Px1-Px3;
       denominator=h3*h2*h1;
       a=(h2*delta3-h3*delta2)/denominator;
       b=(h3^2*delta2-h2^2*delta3)/denominator;
       c=Px3;
%       disp(  '         a                                 b                      c');
%       disp([a,b,c]);
       D=sqrt(b^2-4*a*c);       
       if abs(b-D)<abs(b+D), 
           E=b+D;   		
       else E=b-D;
       end 
       h=-2*Px3/E;
       root=x3+h;
       disp(   '   Iter	                             Aprox-root            Error                             P(root)')
       disp([it,root,abs(h),feval(Poly,root)]);
       if abs(h)==0 | feval(Poly,root)==0
           disp('Final iteration and Exact Solution');
           disp([it,root]);
           return
       end
           if abs(h) <Tol    		
           disp('Final iteration and Approx. root  ');
           disp([it, root]);
           return
       end
       x1=x2;
       x2=x3;
       x3=root;
   end      
   disp('  More iterations are needed'  );