% look at different domains
x1 = linspace(-1, 3, 50);   y1 = linspace(-1, 3, 50);
%x1 = linspace(0, 2, 50);  y1 = linspace(0, 2, 50);
%x1 = linspace(0, 3, 50);  y1 = linspace(0, 2, 50);

[x,y] = meshgrid(x1,y1);

f = (2*x-x.^2).*(2*y-y.^2);

figure(1); clf;
contour(x,y,f,'linewidth',2);
xlabel('x'); ylabel('y');

% optionally, hide f for large values
%f(abs(f) > 3) = nan;

figure(2); clf;
surf(x, y, f);
xlabel('x'); ylabel('y'); zlabel('z');