n <- 2200 X <-runif(n,-1,1); Y <-runif(n,-1,1.5) a <- ifelse( (Y>(X^2)^(1/3)-sqrt(1-X^2))&(Y<sqrt(1-X^2)+(X^2)^(1/3)) ,2,8) x <-runif(n/5,-1,1) ; y <-runif(n/5,-1,1.5) xx <-runif(n/10,-1,1) ; yy <-runif(n/10,-1,1.5) plot(x,y,col=rainbow(50),pch=3) points(Y~X,col=a,pch=3) points(xx,yy,col=rainbow(50),pch=3)