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});