include xpllib; \for Print
def NaN = 1e308;

func real DeCasteljau(C0, C1, C2, T);
real C0, C1, C2, T;
real S, C01, C12;
[S:= 1.0 - T;
C01:= S*C0 + T*C1;
C12:= S*C1 + T*C2;
return S*C01 + T*C12;
];

func real XConvexLeftParabola(T);
real T;
return DeCasteljau(2.0, -8.0, 2.0, T);

func real YConvexRightParabola(T);
real T;
return DeCasteljau(1.0, 2.0, 3.0, T);

func real ImplicitEquation(X, Y);
real X, Y;
return 5.0*X*X + Y - 5.0;

func real F(T);
real T;
real X, Y;
[X:= XConvexLeftParabola(T);
 Y:= YConvexRightParabola(T);
return ImplicitEquation(X, Y) + T;
];

func real Aitken(P0);
real P0;
real P1, P2, P1M0;
[P1:= F(P0);
 P2:= F(P1);
P1M0:= P1 - P0;
return P0 - P1M0 * P1M0 / (P2 - 2.0*P1 + P0);
];

func real SteffensenAitken(PInit, Tol, MaxIter);
real PInit, Tol; int MaxIter;
real P0, P;
int Iter;
[P0:= PInit;
P:= Aitken(P0);
Iter:= 1;
while abs(P-P0) > Tol and Iter < MaxIter do
        [P0:= P;
        P:= Aitken(P0);
        Iter:= Iter+1;
        ];
if abs(P-P0) > Tol then return NaN;
return P;
];

real T0, T, X, Y;
int  I;
[T0:= 0.0;
for I:= 0 to 10 do
        [Print("T0:= %1.1f : ", T0);
        T:= SteffensenAitken(T0, 0.00000001, 1000);
        if T = NaN then
                Print("No answer\n")
        else    [X:= XConvexLeftParabola(T);
                 Y:= YConvexRightParabola(T);
                if abs(ImplicitEquation(X, Y)) <= 0.000001 then
                        Print("Intersection at (%1.6f, %1.6f)\n", X, Y)
                else    Print("Spurious solution\n");
                ];
        T0:= T0 + 0.1;
        ];
]