Here is one approach:
f = @(x,y) x.*y.*(x.^2+y.^2) - 47.3;
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));
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')
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.