The following Mathematica code essentially reproduces the plot:

(* ----- Initialize ----- *)

Clear[x,y,color]

z=x+I y;

xmin=-16;xmax=16;ymin=-16;ymax=16;

W[n_,z_]:=ProductLog[n,z];

color[h_]:={{h,1,h},{h,1,1},{h,h,1},{h,h,h},{1,h,h},{1,0.57+h/2,0},{1,1,h}};

c=color[1/2]; (* light colors *) d=color[0];(* dark colors *)

(* ----- Regions ----- *)

Clear[regn]

regn[n_]:=RegionPlot[Chop[W[n,z Exp[z]]]==z,{x,xmin,xmax},{y,ymin,ymax},PlotStyle->RGBColor[cn+4],PlotPoints->50,FrameLabel->{"x","y"},RotateLabel->False,LabelStyle->(FontSize->10)]

reg=Show[Table[regn[n],{n,-3,3}]];

(* ----- Boundaries ----- *)

Clear[bndn]

bndn[n_]:=ParametricPlot[{-y/Tan[y],y},{y,2ymin,2ymax},RegionFunction->Function[{x,y},(x Cos[y]-y Sin[y])Exp[-x]<=0&& (2n \[Pi]<=y<=(2n+2)\[Pi])],PlotStyle->{RGBColor[dn+4],Thickness[0.015]},PlotRange->{ymin,ymax}]

bnd0=Plot[0,{x,xmin,-1},PlotStyle->{RGBColor[d3],Thickness[0.015]}]; (* xmin < z <= -1 *)

bnd=Show[Table[bndn[n],{n,-3,3}],bnd0];

(* ----- Region Labels ----- *)

Clear[txtn]

txtn[n_,x_,y_]:=Text[Style["n="<>ToString[n],FontSize->15],{x,y}];

xx=-7; tm3=txtn[-3,xx,-15]; tm2=txtn[-2,xx,-10.5]; tm1=txtn[-1,xx,-4]; t0=txtn[0,8,0]; tp1=txtn[1,xx,4]; tp2=txtn[2,xx,10.5]; tp3=txtn[3,xx,15];

txt=Graphics[{tm3,tm2,tm1,t0,tp1,tp2,tp3}];

(* ----- Show all ----- *)

Show[reg,bnd,txt]