Try it this way:
theta = 360 * rand;
rotationMatrix = [cosd(theta) -sind(theta) 0 ;
sind(theta) cosd(theta) 0;
0 0 1];
initialPoints = [ 0 0 0;
15 0 0;
30 0 0;
45 0 0;
60 0 0];
plot3(initialPoints(:, 1), initialPoints(:, 2), initialPoints(:, 3), 'b.', 'MarkerSize', 30);
hold on;
rotatedPoints = initialPoints * rotationMatrix
plot3(rotatedPoints(:, 1), rotatedPoints(:, 2), rotatedPoints(:, 3), 'r.', 'MarkerSize', 30);
grid on;
xlabel('X', 'FontSize', 20);
ylabel('Y', 'FontSize', 20);
zlabel('Z', 'FontSize', 20);
caption = sprintf('The initial points have been rotated by %.2f degrees', theta);
title(caption, 'FontSize', 20);
legend('Initial', 'Rotated');
