Here is one approach:
f = @(x,y) x.*y.*(x.^2+y.^2) - 47.3;
figure
fih = fimplicit(f, [-10 10 -10 10]);
xv = fih.XData;
yv = fih.YData;
v0 = ones(size(xv));
Lpos = xv>=0;
xpos = xv(Lpos);
ypos = yv(Lpos);
xneg = xv(~Lpos);
yneg = yv(~Lpos);
xneg = xneg(~isnan(xneg));
yneg = yneg(~isnan(xneg));
hold on
% fimplicit(f2, [-10 10 -10 10])
yl = ylim;
xl = xlim;
patch([xpos fliplr(xpos)], [ones(size(ypos))*max(yl) fliplr(ypos) ], 'r')
patch([fliplr(xneg) xneg], [ones(size(yneg))*min(yl) (yneg) ], 'r')
hold off
grid
This is not a straightforward problem, so it could be difficult, especially since the NaN values could be difficult to find and elimiinate, and they are not obvious.
The filled areas are those that conform to the stated inequality.
I evaluated them from -10 to 10 with respect to both variables. Change those to get different results.
