> Muller := proc()

local f, x0,x1,x2, a,b,c, dx, D,E, tol;

tol:=10^(-7);

printf("\nEnter a function in terms of variable x:\n");

f := scnf("%a")[1]; f := unaply(f,x);

printf("Enter x0 < x1 < x2:\n");

x0 := scnf("%f")[1]; x1 := scnf("%f")[1]; x2 := scnf("%f")[1];

do

c:=f(x2); # f2

b:=(f(x2)-f(x1))/(x2-x1); # f12

a:=(b-(f(x1)-f(x0))/(x1-x0))/(x2-x0); # f012

b:=b+a*(x2-x1); # (see text)

# Now P(x):=a(x-x2)^2+b(x-x2)+c is the interpolating polynomial

# Rationalize the numerator of {-b +(-) sqrt(b^2-4*a*c)}/(2*a):

D:= sqrt(b^2-4*a*c);

if abs(b+D) > abs(b-D) then E:=b+D; else E:=b-D; fi;

dx := 2*c/E;

x0 := x1; x1 := x2; x2 := x2-dx;

if abs(dx)<tol and abs(f(x2))<tol then

printf(x=x2, ` f(x)` =f(x2)); # x2 may not be real

return();

fi;

od;

end:

> Muller();

Enter a function in terms of variable x:

> x^2+x+1

Enter x0 < x1 < x2:

> 2.5

> 2

> 2.25

x = -.500000000 - .8660254040 I, f(x) = 0