% Muller's Method

f = inline('x^3-5*x^2-6*x-8');

%%% These values find the complex root
%%% -0.59013860341701 + 0.97271617026704i
 x0 = -1;  x1 = 0;  x2 = 1;

%%% These values find the real root 6.18027720683401
% x0 = 4;  x1 = 6;  x2 = 10;

tol      = 10^(-8);
MAX_ITER = 25;

fprintf('---------------------------------\n');
fprintf('  i              x           f(x)\n');
fprintf('---------------------------------\n');

h1 = (x1-x0);
h2 = (x2-x1);
d1 = (f(x1)-f(x0))/h1;
d2 = (f(x2)-f(x1))/h2;
d  = (d2-d1)/(h2+h1);
y0 = f(x0);
y1 = f(x1);
y2 = f(x2);
i  = 3;
while ( i < MAX_ITER )
    b = d2+h2*d;
    D = sqrt(b^2-4*y2*d);
    if abs(b-D)<abs(b+D)
        E = b+D;
    else
        E = b-D;
    end
    h = (-2*y2)/E;
    p = x2+h;
    if (abs(h)<tol)
      fprintf('%3d   % 10.8f   % 12.5e\n',i, x2, y2);
      p0 = p
      return
    else
      fprintf('%3d   % 10.8f   % 12.5e\n',i, x2, y2);  
      x0 = x1;      y0 = y1;
      x1 = x2;      y1 = y2;
      x2 = p;       y2 = f(p);
      h1 = (x1-x0);
      h2 = (x2-x1);
      d1 = (y1-y0)/h1;
      d2 = (y2-y1)/h2;
      d  = (d2-d1)/(h2+h1);
      i  = i+1;
    end
end

% Muller 
%---------------------------------
%  i              x           f(x)
%---------------------------------
%  3    10.00000000    4.32000e+002
%  4    6.15297921   -1.26706e+000
%  5    6.17988202   -1.84868e-002
%  6    6.18027633   -4.11140e-005
%  7    6.18027721    9.49285e-012
% p0 = 6.18027720683401

%---------------------------------
%  i              x           f(x)
%---------------------------------
%  3    1.00000000   -1.80000e+001
%  4   -0.50000000    2.40000e+000
%  5   -0.59938840   -4.75274e-001
%  6   -0.59113078   -5.25936e-003
%  7   -0.59013857    9.10009e-006
%  8   -0.59013860    2.99760e-011
% p0 = -0.59013860341701 + 0.97271617026704i
