I've recently encountered the same problem when trying to implement a Rubik's cube using Golang bindings for RayLib.

Here is my project: https://github.com/geniot/tsp-rubik

My initial approach was the same - use glRotatef to rotate all 3x3x3 cubies around x,y,z axes. As expected my coordinate system got all messed up and cubies were rotating in all possible directions except the ones that I wanted.

Other people on StackOverflow suggested to use quaternions, accumulate rotations into a single something, translate first to origins, then back, use matrices, etc.

But nothing worked for me. It was either too complicated or too vague to implement.

Finally I realized one big truth which I'm going to share with you now :)

It's all about model-view-controller when dealing with OpenGL's transformations. You start with vertices (dots in 3D space) and this is your Model. You use Vertex3f to define your model in the global space.

After that you can start using all fancy ways to manipulate your model. And this is your View - how you want to see your Model. You can transform, scale and rotate but this is not going to change your model!

So back to rotations. When there is a single rotation or the order of rotations is predefined you kind of play a movie to the viewer. This is one thing.

But when you want to give control to the user you should really change the model. At least that's what I find easier in my case.

So our vertices should change their X,Y,Z with each rotation. How do you do that?

I found Vector3RotateByAxisAngle as the ultimate solution for all my rotations. As the comment says in the source code it's based on Euler–Rodrigues formula.

ThreeJS has Vector3#applyAxisAngle.

And this is what I call 'software rotation'. We calculate new positions for our model.

Maybe not the best solution for every use-case but still I hope it helps someone who is facing the same problem.