Mueller's Method (with parabolas)

The following Maple session shows the steps required to find the middle root of the cubic x^3 - 3x +1. The output shown is the result of looping through these steps several times (the loop starts after the initialization of x0, x1, and x2).

> f:=x->x^3-3*x+1;

                                        3
                             f := x -> x  - 3 x + 1
--------------------------------------------------------------------------------
> plot(f(x),x=-2..2);
--------------------------------------------------------------------------------
> x0:=2; x1:=1; x2:=0;

                                     x0 := 2

                                     x1 := 1

                                     x2 := 0
--------------------------------------------------------------------------------
> p:=x->a*x^2+b*x+c;

                                         2
                            p := x -> a x  + b x + c
--------------------------------------------------------------------------------
> eq0:=p(x0)=f(x0);

                                  eq0 := c = 1
--------------------------------------------------------------------------------
> eq1:=p(x1)=f(x1);

             eq1 := .05401353457 a + .2324081207 b + c = .3153288220
--------------------------------------------------------------------------------
> eq2:=p(x2)=f(x2);

             eq2 := .1283921157 a + .3583184557 b + c = -.028950102
--------------------------------------------------------------------------------
> solve({eq0,eq1,eq2});

                   {c = 1., a = .5907265913, b = -3.083276122}
--------------------------------------------------------------------------------
> assign(");
--------------------------------------------------------------------------------
> p(x);

                                    2
                       .5907265913 x  - 3.083276122 x + 1.
--------------------------------------------------------------------------------
> x0:=x1; x1:=x2;

                                x0 := .2324081207

                                x1 := .3583184557
--------------------------------------------------------------------------------
> sol:=evalf(solve(p(x)=0,x));

                         sol := 4.872003013, .347460893
--------------------------------------------------------------------------------
> x2:=sol[2];

                                x2 := .347460893
--------------------------------------------------------------------------------
> a:='a'; b:='b'; c:='c';

                                     a := a

                                     b := b

                                     c := c
--------------------------------------------------------------------------------
>