@acer 

very nice code I wanted exactly this kind of coloring 

I used this article from wikipedia because there is a nice formula for coloring but I have done this only in plot3d like below:

plot3d(abs(x^2 - y^2 + (2*x*y)*I), color = [argument(x^2 - y^2 + (2*x*y)*I), 2*arctan(abs(x^2 - y^2 + (2*x*y)*I))/Pi, 1, colortype = HSV]);

thank you for solving this problem in 2d, do you know maybe what is this colortype function equivalent when we write implicitplot3d?

implicitplot3d(abs(cos(x) + cos(y) + cos(z) + cos(x)*cos(z)*cos(y)*I) = 1, x = -3 .. 3, y = -3 .. 3, z = -3 .. 3, grid = [41, 41, 41], style = surface, color = [sin(x*y), cos(x*y), 0.8, colortype = HSV])

@Kitonum 

I thought about this

https://en.wikipedia.org/wiki/Domain_coloring 

maybe there is some kind of color function in Maple where we can implement code of RGB coloring 

@Carl Love yes, I forgot about it but now it's okay

@acer 

I added two parametrs and I added display to show contours but it doesn't work do you know what is wrong with this modification of your code?

restart;
(m,n) := 200, 200:
(aa,bb,cc,dd) := -3, 3, -3, 3:
M := Array(1..m, 1..n,
           (i,j) -> evalhf( I*(bb-(i-1)*(bb-aa)/(m-1))
                            + (cc+(j-1)*(dd-cc)/(n-1)) ),
           datatype=complex[8], order=Fortran_order):
 

f := (z) -> ((z^2 + a)/(z^3 + b))^3;

P:=proc(MM,param1,param2) local Z,A,temp,img;
  uses ImageTools, plots;
  Z := map[evalhf, inplace](subs(a=param1,b=param2,eval(f)), copy(M)):
  A := Array(1..m, 1..n, 1..3, 1.0, datatype=float[8], order=Fortran_order):
  temp := Array(map[evalhf](z->argument(z), Z), datatype=float[8]):
  A[..,..,1] := map[evalhf,inplace](`*`, FitIntensity(temp,inplace), 360):
  img := HSVtoRGB(A):
  display(plot([[0,0]],x=cc..dd,y=aa..bb,axes=boxed,background=img,
       gridlines=false,size=[400,400]),contourplot(subs(a = param1,b=param2, z = x + y*I, eval(f(z))), x = -3 .. 3, y = -3 .. 3));
end proc:
 

P(M,0.5,2);
 

@Carl Love it's look almost the same (mirror image) but what if function is more complicated you can treat data as variable?

@acer 

restart;
with(plots);
f := z -> z^5 + z^4 - z + 1;
A := plot3d(1, theta = 0 .. 2*Pi, phi = Pi/2 .. Pi, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-f(x + y*I))], axes = none, style = surface, grid = [201, 201]);
B := plot3d(1, phi = 0 .. 2*Pi, theta = 0 .. Pi/2, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-f(I/(y + x*I)))], axes = none, style = surface, grid = [201, 201]);
display(A, B);

g := z -> GAMMA(z - I)^2/GAMMA(z + I)^2;

A1 := plot3d(1, theta = 0 .. 2*Pi, phi = Pi/2 .. Pi, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-g(x + y*I))], axes = none, style = surface, grid = [201, 201]);
B1 := plot3d(1, phi = 0 .. 2*Pi, theta = 0 .. Pi/2, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-g(I/(y + x*I)))], axes = none, style = surface, grid = [201, 201]);
display(A1, B1);
 

@acer 

Thank you, I have done something like this:

S1 := plot3d(1, theta = 0 .. 2*Pi, phi = 0 .. Pi/2, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-g(-I*(y + x*I)))]);
S2 := plot3d(1, theta = 0 .. 2*Pi, phi = Pi/2 .. Pi, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-f(x + y*I))]);
display(S1, S2, orientation = [-10, 57, 52]);
 

and now I see two problems: first is quality of rendering this color function, second is scaling and translate the S2 because rotation is ok and easy to implement we just multiply by some complex number. Do you have any idea how to improve quality and scaling? 

I have also found https://www.wolfram.com/language/12/complex-visualization/riemann-sphere.html?product=language nice post which is about this idea.

@acer 

I used code from https://www.mapleprimes.com/questions/226790-Is-There-Any-Maple-Code-For-The-domain

restart;
with(ImageTools);
m, n := 600, 600;
                        m, n := 600, 600

a, b, c, d := -1.2, 1.2, -1.2, 1.2;
               a, b, c, d := -1.2, 1.2, -1.2, 1.2

M := Array(1 .. m, 1 .. n, (i, j) -> evalhf((b - (i - 1)*(b - a)/(m - 1))*I + c + (j - 1)*(d - c)/(n - 1)), datatype = complex[8], order = Fortran_order);
f := z -> z^5 - z^4 - z + 1;
  f := proc (z) options operator, arrow; z^5-z^4-z+1 end proc

Z := map[evalc, inplace](f, copy(M));
A := Array(1 .. m, 1 .. n, 1 .. 3, 1.0, datatype = float[8], order = Fortran_order);
temp := Array(map[evalhf](z -> argument(-z), Z), datatype = float[8]);
A[() .. (), () .. (), 1] := map[evalhf, inplace](`*`, FitIntensity(temp, inplace), 360);
img := HSVtoRGB(A);
G := plot([[0, 0]], x = c .. d, y = a .. b, axes = boxed, background = img, gridlines = false);
 

mapping 1/z turn inside out this phase portrait so function 1/z^5-1/z^4+1/z+1 give you (after rotation) second phase portrait. 

these plots represent stratographic projection so when you have sphere with upper hemisphere with this phase plot and you project this to plane you get this image but in this situatuion we want to inverse procedure.

@acer