This rotate function works, but results are incorrect because the rotate(X,Y,...) assumes equal scale of X and Y. This is not true for latitude versus longitude. Change in latitude is not the same scale as change in longitude. The rotated object is incorrect. Here is my work around, which works accurately, but requires mapping toolbox functions:
p=geopolyshape(lats,lons);
[distm, azim] = distance(latRef, lonRef, lats, lons, wgs84Ellipsoid);
[latr, lonr] = reckon(latRef, lonRef, distm, azim+angle, wgs84Ellipsoid);
p2 = geopolyshape(latr, lonr);
geoplot(p2), or geoplot(latr, lonr)
Why can't a built-in function do this and make it easy. Or a toolbar, easy access button.