Sorry if the title is a bit cryptic. It's the best I could come up with.
First of all, this question is related to another question I posted here, but that question wasn't posed correctly and ended up generating responses that may be helpful for some people who stumble upon the question, but don't address my problem as it stands.
Also, a disclaimer: I am going to be posting some code here. I will be annotating it to try to make it as clear as possible what I am doing.
Okay, so I am using the Processing programming language to move a box around randomly(according to Perlin noise) As the box moves around, I would like it to rotate such that it looks like it is rolling. I am not looking to simulate exact physics here, I am only looking for a visually pleasing solution.
The general rules for rotation are that, if the box is only moving to the right, then it should rotate around the Y axis such that it rotates to the right and if the box is only moving upwards, it should rotate about the $X$-axis such that it looks like it is rolling directly upwards.
In other words, if this is what the box normally looks like:

Then, if I say: rotate2D(45,0) I should get:

and if I say rotate2D(0,45) I should get:

(for clarity, the function declaration for rotate2D looks like rotate2D(horizontal, vertical) )
Processing makes this easy to do with the rotation functions that it provides, which are:
rotateX();
rotateY();
rotatez();
The problem arises when you want to rotate in 2 dimensions simultaneously. You see, in Processing, If I call rotateX() to rotate an object about the $X$-axis, it is actually the coordinate system itself that gets shifted and all of the other axis go along with it, so when I go to do rotateY() after doing a rotateX()I wont get the result I want because the $Y$-axis that I am rotating around has now been skewed. You can see this in the pictures above.
So, If I do:
rotateX(radians(45));
rotateY(radians(45));
I get:

and if I do:
rotateY(radians(45));
rotateX(radians(45));
I get:

But what I want is something like:

which I drew by calling my function rotate2D()
The problem is that the results that my function produces are a bit counter intuitive. It is not apparent in the above example, but lets say I want to do: rotate2D(90,90);. I would think that would produce a box that looks like It hasn't been rotated at all right? However, it actually produces:

This is my current implementation of rotate2D()
The Function:
void rotate2D(float horizontal, float vertical)
{
float[] rotations ={horizontal, vertical};
if(rotations[0]>rotations[1])
{
float tmp = rotations[1];
rotations[1]=rotations[0];
rotations[0]=tmp;
}
if(rotations[1]==vertical) rotateX(radians(rotations[1]-rotations[0]));
else rotateY(radians(rotations[1]-rotations[0]));
float hIncr= (horizontal<0) ? -0.1 : 0.1;
float vIncr= (vertical<0) ? -0.1 : 0.1;
for(float i=0; i<=abs(rotations[0]);i+=0.1)
{
rotateY(radians(hIncr));
rotateX(radians(vIncr));
}
}
Function Explanation: What the function is doing is essentially as follows (and this is all you really need to address the question from a mathematical standpoint).
The function slowly performs a bunch of $X$ and $Y$ axis rotations right after each other by a small amount each time(0.1 degrees) until it reaches an upper bound. That upper bound is the smaller of the values that the function was passed. The box will have already been rotated about the proper axis by the difference between the larger and the smaller value before this sequence takes place.
To further clarify, I am looking for a mathematical explanation for why this process is producing counter intuitive results, which will hopefully shed light on a solution.
The ideal end goal here is to find a formula that i can use to simulate rotation about a fixed axis in multiple dimensions when all I can do is rotate 'the world' one axis at a time as described above (which is essentially equivalent to rotating the objects relative axes)
please do not hesitate to ask me to clarify further. Thanks in advance.
EDIT
After reading the comments, I now understand that rotations about 2 axis do not commute and I should not expect the results that I was previously expecting.
What I am calculating and passing to my function is the net rotation in each direction.
So, is it possible to edit my formula so that it produces this net rotation? For example, in the 90/90 example, If I continue the sequence for ~37 extra steps (127 total) I get the result that I intuitively want to get) This would need to work for all combinations of inputs though.

Sorry, Luke, I don't think your solution works. You should have a rigorous solution that is founded in the mathematics of rotations. This is most easily done using quaternions--or better yet, their clifford algebra analogues called rotors.
I'll briefly explain clifford algebra rotors and how they can be used to convert sequential rotations into a net rotation.
Clifford algebra
Clifford algebra introduces a product of vectors that is noncommutative but still associative. If you're familiar with dot and cross products, it incorporates properties of both. Given an orthogonal basis $e_1, e_2, e_3$, we have the following:
$$e_i e_j = \begin{cases} +1 & i = j \\ -e_j e_i & i \neq j \end{cases}$$
Again, it's also associative, so $e_1 e_2 e_3 = (e_1 e_2) e_3 = e_1 (e_2 e_3)$ for instance. This creates an 8-dimensional vector space, with the following basis elements: $$1 \\ e_1, e_2, e_3 \\ e_1 e_2, e_2 e_3, e_3e_1 \\e_1 e_2 e_3$$
Scalar multiples of $1$ are scalars. Linear combinations of $e_1, e_2, e_3$ are vectors. Linear combinations of $e_1 e_2, e_2 e_3, e_3 e_1$ are called bivectors, and these are of interest for rotations. Rotations take place in planes (it's only in 3d that we can say they are about an axis instead, equivalently), and bivectors directly describe the planes in which rotations take place.
Rotors and rotations
I'll state the following without proof, though it should be familiar to anyone with a base understanding of quaternions.
Define a rotor as the exponential of a bivector (usually defined through a power series). Let $B$ be some unit bivector. Then the rotor $q = \exp(Bt)$ for some scalar $t$ is
$$q = \exp(Bt) = \cos t+ B \sin t$$
A rotation map $\underline R$ in the plane corresponding to $B$ by the angle $\theta$ takes the form
$$\underline R(a) =e^{-B\theta/2} a e^{B\theta/2}$$
where $a$ is any vector.
(Note: for those with quaternion background, $i = -e_2 e_3$, $j = -e_3 e_1$, and $k = -e_1 e_2$, so there is no sign discrepancy with this definition.)
Application: finding a net rotation
Using rotors simplifies a lot of the mathematics of analyzing rotations. For instance, let's consider the problem you have: your rotation functions rotate the global frame, and themselves use the global frame as the frame of reference for these rotations. This makes sequencing the rotations somewhat nontrivial, but rotors let us find a simpler way of doing things.
Let $e_x, e_y, e_z$ be the original frame. If you want to rotate about the original $e_x$ and then about the original $e_y$, then your rotors would look like this:
$$q_\text{net} = \exp(-e_ze_x \phi/2) \exp(-e_y e_z \theta/2) = q_y q_x$$
We can then multiply this out to find the net rotation plane and the net angle.
$$q_\text{net} = \cos \frac{\phi}{2} \cos \frac{\theta}{2} - e_z e_x \sin \frac{\phi}{2} \cos \frac{\theta}{2} - e_y e_z \cos \frac{\phi}{2} \sin \frac{\theta}{2} + e_x e_y \sin \frac{\theta}{2} \sin \frac{\phi}{2}$$
We can find the net rotation angle $\alpha$ by
$$\cos\frac{\alpha}{2} = \cos \frac{\phi}{2} \cos \frac{\theta}{2}$$
Let's consider the case like you suggested, $\phi = \theta = \pi/2$. Then we should have
$$\cos\frac{\alpha}{2} =\left (\cos\frac{\pi}{4}\right)^2 = \frac{1}{2} \implies \frac{\alpha}{2} = \frac{\pi}{3}$$
or $\alpha = 2\pi/3 = 120^\circ$.
The net rotation plane is $(e_y e_z + e_z e_x - e_x e_y)/\sqrt{3}$, or a rotation about the axis $(e_x + e_y - e_z)/\sqrt{3}$.
Edit: with this in mind, here's what you need to do.