function result = muller ( f, x0, x1, x2, FATOL, XRTOL, ITMAX)
%
%  MULLER carries out Muller's method for seeking a real root of a nonlinear function.
%  This is a "stripped down" version with little error checking.
%
% original written by John Burkardt
% 
% Thu Oct 13 23:40:41 EDT 2005
% Modified by Andy Long, to comply with Newton's form of the interpolating
% polynomial, using divided differences; also to compute complex roots. Note
% that this algorithm works with arbitrary functions!
%
y0 = feval ( f, x0);
y1 = feval ( f, x1);
y2 = feval ( f, x2);

% Check to see if we've got a root at the outset:
if ( y0 == 0 )
  result = x0;
  return
elseif ( y1 == 0 )
  result = x1;
  return
elseif ( y2 == 0 )
  result = x2;
  return
end

% Compute some differences:
x10=(x1-x0);
x21=(x2-x1);
y10=(y1-y0);
y21=(y2-y1);
% divided differences:
f01 = y10/x10;
f12 = y21/x21;
f012 = (f12-f01)/(x2-x0);

it = 0;
format long

while ( it <= ITMAX ) 
  
  it = it + 1;

  % Determine the coefficients A, B, C 
  % of (Newton's form of) the polynomial 
  % which goes through the data (X0,Y0), (X1,Y1), (X2,Y2).

  % Newton's form:
  % Y(X) = A*(X-X2)(X-X1) + B*(X-X2) + C
  % add the appropriate form of zero:
  % Y(X) = A*(X-X2)(X-X2+X2-X1) + B*(X-X2) + C
  % regroup:
  % Y(X) = A*(X-X2)(X-X2) + (B+A*(X2-X1))*(X-X2) + C

  a = f012;
  b = f12+a*x21;
  c = y2;

  if ( a ~= 0 ) % we're really quadratic:
    
    disc = b * b - 4.0 * a * c;
    
    % Here we check to see which root is closer to the approximate root x2:
    % choose the roots of the quadratic polynomial, and choose the one that makes
    %the LEAST change to X2
    q1 = ( b + sqrt ( disc ) );
    q2 = ( b - sqrt ( disc ) );
    if ( abs ( q1 ) < abs ( q2 ) )
      dx = - 2.0 * c / q2;
    else
      dx = - 2.0 * c / q1;
    end
    
  elseif ( b ~= 0 ) % We're linear:
    dx = - c / b;

  else % We're constant! Whooooooooooooooooooooaaaaaaaaaa.....
    '(Muller algorithm broke down, results unreliable.)'
    result = x2;
    return
  end
  
  tmp = x2 + dx;
  x0 = x1;
  x1 = x2;
  x2 = tmp;  
  y0 = y1;
  y1 = y2;
  y2 = feval ( f, x2);  

  % Report the current iteration:
  [ x2, y2 ]

  % Declare victory if the most recent change in X is small, and 
  % the size of the function is small:
  if ( abs ( dx ) < XRTOL * ( abs ( x2 ) + 1.0 ) & abs ( y2 ) < FATOL ) 
    result = x2;
    return
  end

  % Otherwise, update! Move values down in order where possible, and
  % calculate the new:
  x10=x21;
  x21=(x2-x1);
  y10=y21;
  y21=(y2-y1);
  f01 = f12;
  f12= y21/x21;
  f012 = (f12-f01)/(x2-x0);
end

% If maxits has been exceeded, then 
result = x2;
'(Maximum number of iterations taken, results unreliable.)'
endfunction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
'Some examples: notice that only one is a polynomial.... ael'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
'zero of sine:'
muller("sin",.1,.2,.3,1e-6,1e-6,10)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
'zero of cosine:'
muller("cos",.1,.2,.3,1e-6,1e-6,10)
'true value: pi/2:'
pi/2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
'finding real and complex roots: (x-1)*(x^2+1)'
function f = f(x)
f=(x-1)*(x^2+1);
endfunction
'Looking for that real root...'
muller("f",.1,.2,.3,1e-6,1e-6,10)
'Looking for a complex root...'
muller("f",.1,.2,.3+i,1e-6,1e-6,10)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
'Another non-polynomial function: exp(x)-21'
function f = f(x)
f=exp(x)-21;
endfunction
muller("f",.1,.2,.3,1e-17,1e-17,10)
feval("log",21.0)