That's what happens when you run contour() on data with step discontinuities. If the step goes from Zmax to Zmin, then it will contain all contour lines, nonintersecting.
Since you're not actually using the contours, you might be able to just use pcolor() or image()/imagesc().
Consider the example:
x = 1:100;
y = x';
z = x+y;
z = circshift(z,[0 -30]); % create step discontinuities
nlevels = 12;
% plot using contourf
[~,hc] = contourf(x,y,z,nlevels-1);
hc.LineStyle = 'none';
colormap(hsv(nlevels))
Note the behavior at the step discontinuities
% plot using pcolor instead
figure
pcolor(x,y,z)
colormap(hsv(nlevels))
shading flat
The disadvantage to doing things this way is the lack of edge interpolation. If the data resolution is low, it might be quite pronounced.
EDIT:
I should note the banding at step discontinuities is reduced somewhat for larger arrays; however, attempting to upscale existing data with any continuous interpolant will spread out the discontinuity by at least the same factor the data is upscaled. You either have to start with fine data or use nearest-neighbor interpolation (which will make the contours jagged like using pcolor()).
