f1=x^3 + y^3 - 2 x y;
f2=x^2 + y^3-1;

{r1,r2,r3}={x,y}/.Solve[{f1==0,f2==0},{x,y},Reals]//N

jac=D[{f1,f2},{{x,y}}];

newton[{x0_,y0_}]:=({x,y}-Inverse[jac].{f1,f2})/.{x->x0,y->y0}

barva[{x0_,y0_}]:=(n=1;xy={x0,y0};e=10^-9;While[
xy=newton[xy];
d1=Norm[xy-r1];d2=Norm[xy-r2];d3=Norm[xy-r3];
(Min[d1,d2,d3]≥e)&& (n++)<20];
b=Black;
If[d1<e,b=Red];
If[d2<e,b=Green];
If[d3<e,b=Blue];
b);

barva[{1.3,1.5}]

ArrayPlot[Table[barva[{xw,yw}],{xw,-3.0021,3,0.01},{yw,-3.0062,3,0.01}]]