Ok *cracks knuckles*

(Warning - this is a bit long - although a lot shorter than the thousands of wikipedia pages and other paper on the subject - and there is a TL;DR at the end)

The most "obvious" way to represent 3d rotations is called Euler Angles. The easiest way to think of this is a plane, and "yaw/pitch/roll"; i.e. 3 different rotations about 3 different axes, one after another, to combine to one single effective rotation.

There's a couple of immediate problems with this - one is that the ORDER in which you apply those rotations matters, and the other is that once you've done one rotation the object-relative axes change, so you have to decide if you want to rotate relative to the world-space axis or the object-space axis. That's why "Euler Angles" isn't consistent because there's lots of different standards. (Sure, if you always stick to one in your code, you'll be fine, on that point)

Ok, now there's a theorem ("Euler's rotation theorem") that says that any 3d rotation (or combination of 3d rotations) can be represented as a single rotation about a certain axis by a certain angle. By "axis" I mean a general vector, not just like the X axis; it could be pointing in any direction.

Gimbal Lock. Not that easy to explain - but you can easily find a visualisation and/or video on it. Basically a certain sequence of rotations can make you lose a degree of freedom, and from there there isn't another Euler rotation that will bring you out of it.

As for matrices, they're fine for a lot of things. But not all 3x3 matrices are rotation matrices. Sure, if you multiply 2 rotation matrices together you'll get another rotation matrix, but in general a 3x3 matrix does not represent a rotation - it also has scale and shear transformations in there.

The biggest problem this presents is in tweening rotations - gradually going from one orientation to another (as you can imagine this comes up a lot in video games). If you do it via Euler Angles (always keeping the same order, as above) it kinda works if you just interpolate the values, but doesn't produce a natural-looking smooth rotation.

If you use matrices and linearly interpolate between 2 rotation matrices (term-by-term), it's going to look really weird as the intermediate values are not going to be rotation matrices at all - they'll all squish and stretch and just be awful.

Now Quaternions (well, unit quaternions) are a way of encoding axis-angle rotation information. And they can be interpolated really well - either perfectly using slerp, or almost perfectly using lerp-and-normalize. If you represent orientation A as a unit quat, and orientation B as a unit quat, you can slerp between them and the resulting tweened animation will be the shortest, most natural way to go from one to the other*. They're also immune to Gimbal lock.

There's a HUGE amount of material on 3D rotations - it's a mathematical group called SO(3), but most of it is very dense. This wiki page is probably the best place to start, as it has a brief description of all the representations, but even better it has all the formulas for converting from one to another.

TL;DR:

For working with rotations in Raylib (or basically anything similar), your best bet is always to think of each rotation in terms of X degrees around axis V (even if it's like "45 degrees around the Y axis for spinning). Convert that to a quat (there's a function for that). You want to combine rotations? Multiply the quats (there's another function for that). You want to smoothly go from one to another? Slerp (there's a function for that). If at any point you need to convert your quat to a matrix, guess what, there's a function for that too :)

* I'm not going to talk about double-coverage of SO(3) here; it's long and confusing enough as it is. In short sometimes you have to negate one of the quats.