Метод Мюллера в Mathematica

Question

Я хотел бы написать функцию, в которой есть цикл, который выполняет операции, необходимые для метода Мюллера .

f[x_] := x^3 - x - 1;
x0 = 0.8
x1 = 1.5
x2 = 2.0
x3 = 5.0;
\[Epsilon] = 0.001;

While[(Abs[f[x3]] >= \[Epsilon]),
 h0 = x1 - x0;
 h1 = x2 - x1;
 d0 = (f[x1] - f[x0])/h0;
 d1 = (f[x2] - f[x1])/h1;
 A = (d1 - d0)/(h1 + h0);
 B = A*h1 + d1;
 Cx = f[x2];
 raiz = Sqrt[B^2 - 4.0*A*Cx];
 If[Abs[B + raiz] > Abs[B - raiz], dens = B + raiz, dens = B - raiz];
 x3 = (x2 - 2*Cx)/dens;
 i++;
 Print["Iteration: ", i, "\t root \[TildeTilde] ", x3];
 x0 = x1;
 x1 = x2;
 x2 = x3;
 ]

Но я получаю бесконечный цикл ...

Answer 1

метод Мюллера по Эрику (всегда лучше, чем в Википедии): Спасибо Хайке за указание на несколько ошибок в комментарии ниже

h[x_] := HermiteH[24, x];
i = Length@CoefficientList[h[x], x] - 1;
f[i, x_] := h[x];
roots = {};
While[ i > 1,
  x0 = -2; x1 = -1; x2 = -.5; k = 1;
  While[Abs[k] > .001,
   q = (x0 - x1)/(x1 - x2);
   a = q f[i, x0] - q (1 + q) f[i, x1] + q^2 f[i, x2];
   b = (2 q + 1) f[i, x0] - (1 + q)^2 f[i, x1] + q^2 f[i, x2];
   c = (1 + q) f[i, x0];
   p = Sqrt[b b - 4 a c];
   xp = x0 - (x0 - x1) 2 c /(k = If[Abs[b + p] > Abs[b - p], b + p, b - p]);
   {x2, x1, x0} = {x1, x0, xp};
   ];
  AppendTo[roots, xp];
  i--;
  f[i, x_] = f[i + 1, x]/(x - xp);
  ];
Show[
 Plot[h[x], {x, -2, 2}],
 Graphics[{PointSize[Large], Point[{#, 0} & /@ roots]}]]