This is how I would do it:
F = @(x,y) [1*(x>=1 & y>=1) + x.*((x>=0 & x<1) & (y>=1)) + y.*((x>=1) & (y>=0 & y<1)) + x.*y.*((x>=0 & x<1) & (0<=y & y<1))];
[X,Y] = meshgrid(linspace(-1, 3, 50));
FXY = F(X,Y);
figure(1)
surfc(X, Y, FXY)
grid on
xlabel('X')
ylabel('Y')
zlabel('F(X,Y)')
view([130 40])
producing this plot:
