@lekcha I have a code 

restart: 

newton := proc(x, y)
    local z, m;
    z := evalf(x+y*I);
    for m from 0 to 15 while abs(z^3-1) >= 0.001 do
       z := z - (z^3-1)/(3*z^2)
    od;
    m
 end:

P := plots:-densityplot(newton, -2 .. 2, -2 .. 2, grid = [500, 500], colorstyle = HUE, axes = none, scaling = constrained, style = patchnogrid, size = [500, 500])
applyop(`*`, [1,4,2], P, .85);

got the picture as follows

If I want to assign a color to an image based on the m variable in the Newton function, for example:

cm := ["turquoise","cyan","blue", "brown","coral", "maroon", "orange", "pink", "magenta", "red", "sienna", "violet", "wheat", "yellow", "aquamarine", "black"]:

or any color but m=15 items are black

how should i fix it

@acer Thank you, I've done it.

@acer I'm sorry that my English misunderstood you.

 i have

m:= Matrix([[1, 2, 3], [4, 5, 6]]):

cm := ["aquamarine","blue","brown",
        "coral","cyan","black"]:

to show

this is what i want.

@acer thanks for the answer.

I'm going to make an image similar to Fractals, but the difference is that each point in the image is the number of convergence cycles. which I came up with a new method thinking that it can be created as a matrix where members are converging cycles If I had that matrix How can I create an image?
m:= Matrix([[1, 2, 3], [4, 5, 6]]):

cm := ["aquamarine","blue","brown",
        "coral","cyan","black"]:

Sorry my English is not good.

@acer what you offer It's exactly what I want.
But when I follow the code It gives an error "Error, HeatMap is not a command in the Statistics package", what should I do? will get the picture as you do

@acer I tried it and got an error. Error, HeatMap is not a command in the Statistics package. 

How should I solve the problem?

maple 18

@acer thank you very much you made me gain more knowledge

COR:= proc(f, x0, tol)
local y, z, s, i, fs, a, de, n;
s := x0:
a[0]:=x0;
 fs := f(x0);
  de := 2*tol;
for i from 1 to 10^6
while abs(fs)>tol and de>tol do
y := evalf(s- f(s)/f1(s)):
z:=evalf(s-(f(s)+f(y))/f1(s)):
s:=evalf(z-(f(y)^2*(2*f(s)+f(y)))/(f(s)^2*f1(s))):
fs:=evalf(f(s));
a[i]:=s:
de:=abs(evalf(a[i]-a[i-1]));
n:=i:
end do:
return n, s, fs, de;
end proc;

f := x -> x*exp(x^2)-sin(x)^2+3*cos(x)+5:
f1 := D(f):
f2:=D(f1):
x0:=-1:
Digits :=100: 
tol := 10.^(-16):

sol := COR(f, x0, tol);
%;
Error, (in COR) cannot determine if this expression is true or false: 0. < -exp(1)-sin(1)^2+3*cos(1)+4.999999999999999900000000000000000000000000000000000000000000000000000000000000000000000000000000000
 

I tried another equation and it worked pretty well. But this equation shows an error. Where do I need to fix it?

@acer Thank you, it helps a lot.