Help:=proc(): print(``):end:


f := x -> x^3 - 44*x^2 + 564*x-1728; 
Digits := 15; 
epsilon := 10^(-6);
Df := D(f);

na:= x -> evalf(x-f(x)/Df(x));

nit := proc(xinit)
local xcurrent,xnext;
	xcurrent := xinit;
	xnext := na(xinit);
	while abs(f(xcurrent)-f(xnext)) > epsilon do
		xcurrent := xnext;
		xnext := na(xnext);
	od;
xnext;
end proc;

nit(1); nit(14); nit(30);

solve(f(x));
rho1,rho2,rho3 := %;

simplify(rho1+rho2+rho3);
simplify(rho1*rho2+rho2*rho3+rho3*rho1);
simplify(rho1*rho2*rho3);

simplify((1/rho1^2)+(1/rho2^2)+(1/rho3^2));

simplify((1/rho1^2)+(1/rho2^2)+(1/rho3^2),
[r1+r2+r3=e1,r1*r2+r2*r3+r3*r1=e2,r1*r2*r3=e3]);

rho1,rho2,rho3 := fsolve(f(x));
BArho1 := {}:
BArho2 := {}:
BArho3 := {}:

square := 30;
step := 0.5;

for xi from -square by step to square do
for yi from -square by step to square do

	x0 := evalf(xi+yi*I);
	aux := nit(x0);

if abs(aux - rho1) < epsilon then BArho1 := BArho1 union {x0}; end if;
if abs(aux - rho2) < epsilon then BArho2 := BArho2 union {x0}; end if;
if abs(aux - rho3) < epsilon then BArho3 := BArho3 union {x0}; end if;

od;
od;

nops(BArho1);nops(BArho2);nops(BArho3);

%+%%+%%%;    (2*square*(1/step)+1)^2;

with(plots):

BA1Plot := complexplot(BArho1,style=point,axes=box,color=red):
BA2Plot := complexplot(BArho2,style=point,axes=box,color=green):
BA3Plot := complexplot(BArho3,style=point,axes=box,color=blue):

display({BA1Plot,BA2Plot,BA3Plot});