I can see several ways off the top of my head.
Consider some arbitrary implicit function.
fun = @(x,y) x.^2 + y.^3 - 3;
Easy enough to plot the region boundary, so the locus of where it equals zero.
But which half of that domain is the negative part? pcolor should work.
[X,Y] = meshgrid(linspace(-4,4,500),linspace(-3,2,500));
Fxy = fun(X,Y);
H = pcolor(X,Y,double(Fxy < 0);
H.EdgeColor = 'none';
colormap([1 1 1;.1 .6 .7])