function x=muller(FUN, x0,x1,x2,TOL,N);
%%
%%   function x=muller(FUN, x0,x1,x2,TOL,N);
%%   to find a solution to f(x)=0 given three approximations x0,x1,x2;
%%
%%   INPUT:
%%   x0,x1,x2: the initial values;
%%   TOL: tolerence;
%%   N: maximum number of iteration
%%   FUN: the given function;
%%
%%   OUTPUT:
%%   approximation solution x or massage of failure;
h1=x1-x0;
h2=x2-x1;
y0=feval(FUN,x0);
y1=feval(FUN,x1);
y2=feval(FUN,x2);
g1=(y1-y0)/h1;
g2=(y2-y1)/h2;
d=(g2-g1)/(h1+h2);
i=2;
while i<=N
   i
   b=g1+h2*d;
   u=sqrt(b.^2-4*y2*d);
   if abs(b-u)<abs(b+u)
      v=b+u;
   else
      v=b-u;
   end;
   h=-2*y2/v;
   x=x2+h
   if abs(h)<TOL
      break;
   end;
   x0=x1;
   x1=x2;
   x2=x;
   y0=feval(FUN,x0);
   y1=feval(FUN,x1);
   y2=feval(FUN,x2);  
   h1=x1-x0;
   h2=x2-x1;
   g1=(y1-y0)/h1;
   g2=(y2-y1)/h2;
   d=(g2-g1)/(h1+h2);
   i=i+1;
end;
if i>N
   massage='method failed aftern0 iteration'
end;      
