You have a surface, generated using meshgrid. I would split that surface into a triangulated mesh. So two triangles per rect in the mesh. Next, for each triangle in the triangulation, I would now generate a normal vector to the surface, using the function cross. Test which direction that normal vector points, up or down. Essentially, you will just use the z-component of the surface normal. Color the surface depending on the direction identified by the normal vector. Simple.
For example, a simple paraboloid of revolution.
[t,x] = meshgrid(linspace(0,2*pi),linspace(0,10,50));
y = sqrt(x).*cos(t);
z = sqrt(x).*sin(t);
% triangulate the grid
[m,n] = size(t);
[in,im] = meshgrid(1:n-1,1:m-1);
ind = sub2ind([m,n],im(:),in(:));
tri = [[ind,ind+1,ind+m];[ind+1,ind+m+1,ind+m]];
xyz = [x(:),y(:),z(:)];
normals = cross(xyz(tri(:,2),:) - xyz(tri(:,1),:),xyz(tri(:,3),:) - xyz(tri(:,1),:));
direction = normals(:,3) >= 0;
H1 = trimesh(tri(direction,:),xyz(:,1),xyz(:,2),xyz(:,3));
hold on
H2 = trimesh(tri(~direction,:),xyz(:,1),xyz(:,2),xyz(:,3));
hold off
H1.FaceColor = 'r';
H2.FaceColor = 'b';
H1.EdgeColor = 'none';
H2.EdgeColor = 'none';
I chose a coloring that has faces in the upper branch as blue, those on the lower branch are red.
Actually finding the fold in the curve itself seems, at least initially, a more difficult problem. You could view that as an implicit function to locate the path of a singularity in the surface. As I think though, perhaps simplest is to locate triangles in the triangulation I generated where the triangles share an edge, but the normal vectors point in different directions. The shared edge will delineate an approximate path around the surface of the "fold". As such, it would be quite doable.
