This is an interesting problem!
This should get you started —
numPatches = fix(Facets/5);
[X,Y,Z] = sphere(Facets-1);
hs = surf(X, Y, Z, 'EdgeColor','none');
idxk = 1:fix(Facets/numPatches);
idx = idxk + numel(idxk)*(k-1);
hs.CData(:,idx) = k-fix(Facets/2);
NOTE — The surface patches are defined by the CData colormap colors that go from -1 to +1 so whatever the desired colormap is (I chose turbo that was introduced in R2021a) will automatically map to those values. (After a false start using patch, this turned out to be reasonably straightforward.) The loop works by selecting the columns defined by ‘idx’ and then colouring them according to the value of the colormap at that point. The rotate call is optional, and simply inclines the sphere object so that more of its patches are visible.
Experiment to get different results.