Draw Coordinate System with patch() command

Hey folks,

i need to visualize the spacial rotation of a cuboid efficiently. For this task, i use the patch() command as shown in the slides of the package: Using patch() and rotate() - Basics, which is probably the way to go for me, since it allows to draw frames efficiently. Im already able to draw and rotate the cuboid. To make my visualization easier to interpret, I would like to draw two orthogonal coordinate systems now. The first coordinate system represents the inertial coordinate system. The second coordinate system is fixed to the cuboid and rotates with it. Can someone help me and tell me how to draw the coordinate systems? If possible, I would like to label the coordinate system axes or at least draw them in different colors.

Here is a minimum example of my code:

% define array of edge point coordinates [x; y; z]
edgePoints = [..];

% define six surfaces of cuboid 
surfaces = [..];

% create figure
figure();

% create patch
p = patch('faces',surfaces,...
          'vertices',edgePoints,...
          'facecolor',[.5 .5 .5],...
          'edgecolor',[1,1,1],...
          'facealpha',0.6);

% in this loop, the cuboid gets rotated     
for i = 1 : numberOfTimeSteps
    
    V = get(p, 'vertices');
    
    % rotate edge pints
    for j = 1 : size( edgePoints, 1)
        V(j,:) = ( RotationMatrix(:,:,i)*edgePoints(j,:)' )';
    end
    
    set(p,'vertices', V);
    drawnow
   
end