function [k,x,y,err,S,F]=muller(f,x2,x0,x1,delta,epsilon,max1)
%Input - f is the object function input as a string 'f'
%	x0, x1, and x2 are the initial approximations 
%      - delta is the tolerance for x0, x1, and x2
%      - epsilon the the tolerance for the function values f 
%      - max1 is the maximum number of iterations
%Output - k is the number of iterations that were carried out
%       - x is the Muller approximation to the zero of f
%       - y is the function value y = f(x)
%       - err is the error in the approximation of x.
%       - S' contains the sequence {x}
%       - F' contains the sequence {f(x)}

%  NUMERICAL METHODS: Matlab Programs
% (c) 2004 by John H. Mathews and Kurtis D. Fink
%  Complementary Software to accompany the textbook:
%  NUMERICAL METHODS: Using Matlab, Fourth Edition
%  ISBN: 0-13-065248-2
%  Prentice-Hall Pub. Inc.
%  One Lake Street
%  Upper Saddle River, NJ 07458

format short;
%format long;
disp('iteration      x2         x0        x1      f(x0)')
%Initalize the matrices X and Y
X=[x2 x0 x1];
y=feval(f,x0);
D=[0,X,y];
disp(D); 
Y=feval(f,X);
for k=1:max1
   h1=x1-x0;
   h2=x0-x2;
   G=h2/h1;
   c=Y(2);
   a=(G*Y(3)-Y(2)*(1+G)+Y(1))/(G*h1*h1*(1+G));
   b=(Y(3)-Y(2)-a*h1*h1)/h1;
	%Suppress any complex roots
   	if b^2-4*a*c > 0
      	disc=sqrt(b^2-4*a*c);
   	else
      	disc=0;
   	end
   	%Find the closest  root 
   	if b < 0
      	disc=-disc;
   	end
   z=2*c/(b+disc);
   x=x0-z;
   	if x > x0
      	x2=x0;
      	x0=x;
   	else
      	x1=x0;
      	x0=x;
   	end	
   S(k)=x;
   X=[x2 x0 x1];
   Y=feval(f,X);
   y=feval(f,x);
   F(k)=y;
   D=[k,X,y];
   disp(D); 
   %Determine stopping criteria
   err=abs(z);
   relerr=err/(abs(x)+delta);
   if (err<delta)|(relerr<delta)|(abs(y)<epsilon)
      break
   end
end
S=S';
F=F';
