Rotate an object on its own axes in OpenGL

So I have a cube that I wish to rotate around any of its three axes (the axes of the cube, not the window). As many other similar questions have stated, my rotations work as long as I am only rotating in one direction, but when I start mixing them, I get strange results. In particular, the rotation about the y-axis always rotates about the y-axis of the window, regardless of how the cube has been rotated.

My drawing code is as follows:

glMatrixMode(GL_MODELVIEW) ;
glLoadIdentity();
gluLookAt(0.0, 5.0, 15.0,  
          0.0, 0.0, 0.0,  
          0.0, 1.0, 0.0);
glPushMatrix();
glRotatef(theta_y,0,1,0); 
glRotatef(theta_x,1,0,0);
glRotatef(theta_z,0,0,1);
draw_cube();
glPopMatrix();

This question seems to describe pretty much exactly what I am trying to do, and the accepted answer seems to be what I want to do, however the link he provides is dead.

From what I can gather in the linked question, my z-rotation is performed first, which means in my x-rotation (which would be next) I need to, instead of rotating about (1,0,0), I would rotate about (-sin(theta_z),cos(theta_z),0). But then for the third rotation, he only gives the link, saying it gets "very complicated." Now that the link is dead, I'm not sure how to go about this third rotation.

EDIT: Based on the replies, I have added three vectors: cube_x, cube_y, and cube_z to hold the current axes of my cube.They are initialized as follows:

float cube_x[] = {1.0f,0.0f,0.0f,0.0f};
float cube_y[] = {0.0f,1.0f,0.0f,0.0f};
float cube_z[] = {0.0f,0.0f,1.0f,0.0f};

Back in my drawing function, I have changed it so that the rotations are done around these axes rather than the global axes as I had previously. After performing the three rotations, I call glGetFloatv to get the current matrix and use a new function, update_vector to fill my cube vectors with their new values. This all is included below:

void update_vector(float vector[],float matrix[]) {
    vector[0] = matrix[0]*vector[0] + matrix[1]*vector[1] + matrix[2]*vector[2] + matrix[3]*vector[3];
    vector[1] = matrix[4]*vector[0] + matrix[5]*vector[1] + matrix[6]*vector[2] + matrix[7]*vector[3];
    vector[2] = matrix[8]*vector[0] + matrix[9]*vector[1] + matrix[10]*vector[2] + matrix[11]*vector[3];
    vector[3] = matrix[12]*vector[0] + matrix[13]*vector[1] + matrix[14]*vector[2] + matrix[15]*vector[3];

}

void my_display(void) {
    glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT) ;

    glMatrixMode(GL_MODELVIEW) ;
    glLoadIdentity();
    gluLookAt(0.0, 5.0, 15.0,  
          0.0, 0.0, 0.0,  
          0.0, 1.0, 0.0); 
    glPushMatrix();
    glRotatef(theta_y,cube_y[0],cube_y[1],cube_y[2]); 
    glRotatef(theta_x,cube_x[0],cube_x[1],cube_x[2]);
    glRotatef(theta_z,cube_z[0],cube_z[1],cube_z[2]);

    //get the current matrix
    float my_matrix[16];
    glGetFloatv(GL_MODELVIEW_MATRIX, my_matrix);

    //Multiply the matrix by each of my vectors
    update_vector(cube_x,my_matrix);
    update_vector(cube_y,my_matrix);
    update_vector(cube_z, my_matrix);

    make_cube();
    glPopMatrix();

    /* buffer is ready */
    glutSwapBuffers();

    return ;
}

Now when I try to rotate about the y axis, the cube seems to wobble around along a few axes for a little bit before finally settling into a rotation about the x-axis. Rotating about Z does the same thing. Only rotating about the x-axis seems to work.

EDIT: SOLUTION:

In his comment below, Gavin mentioned that the Y-axis rotation was always happening first. This is the whole issue. By adding a flag ROT and setting it based on which axis I am attempting to rotate around, I am able to order the rotations so that the axis being rotated about is done last.