To rotate a vector around an axis by an angle 'a', you multiply the vector with the appropriate matrix:

Rotate around X axis:

[ 1 0 0 0 ]

[ 0 cos a sin a 0 ]

[ 0 -sin a cos a 0 ]

[ 0 0 0 1 ]

Rotate around Y axis:

[ cos a 0 -sin a 0 ]

[ 0 1 0 0 ]

[ sin a 0 cos a 0 ]

[ 0 0 0 1 ]

Rotate around Z axis:

[ cos a sin a 0 0 ]

[-sin a cos a 0 0 ]

[ 0 0 1 0 ]

[ 0 0 0 1 ]

You need to watch out for what is known as "Gimbal lock", where one or more of the rotations have no effect due to the order in which you perform the rotation.

To find the angle between two vectors around each axis, treat each vector as being 3 different 2D vectors, and use the dot product:

For vectors u and v:

u dot v = u.x * v.x + u.y*v.y + .. = |u| * |v| * Cos a

Therefore:

a = arccos( u dot v / (|u| * |v|) )