I see, if no boundary, it will be easy with the help of built-in functions:
1. [x,y] = meshgrid(input arguments)% get the mesh grid point
2. z = (x-3).^2 + (y-3).^2;%calculate the z value
3. surface(x,y,z);% plot the surface
4. shading interp% make the surface more pretty
I get quite confused about your boundary condition: g = 10(x1+x2≥4); maybe you should set z value at these points to NaN, and continue with surface(x,y,z); shading interp
z(points out the boundary) = nan;
