Assume you have a vector (x,y,z) and you want to rotate it to, say x-axis, you can multiply the rotation matrix to the vector.
First, make the vector a column vector and append 1 to the end. It becomes a 4×1 matrix:
v = x y z 1
Then, convert the vector from Cartesian to spherical coordinate:
[theta, phi, r] = cart2sph(x,y,z)
The meaning of theta, phi and r is illustrated by the following diagram:

Then based on theta and phi, do proper rotations. For example, to rotate the vector to x-axis, you might want to rotate the vector to XZ plane first (by rotate along Z axis with angle theta), and then to x-axis by rotate along Y axis with angle phi. You need to construct the rotation matrix (below) and multiply the matrix with the original vector. Note the order matters.
v_final = M2*M1*v
The following material is from http://www.siggraph.org/education/materials/HyperGraph/modeling/mod_tran/3drota.htm
Rotate along z-axis

         ( cos q  sin q  0  0)
         (-sin q  cos q  0  0)
         ( 0        0    1  0)
         ( 0        0    0  1)
Along x-axis

        (1    0      0    0)
        (0  cos q  sin q  0)
        (0 -sin q  cos q  0)
        (0    0     0     1)
Along y-axis

        (cos q  0  -sin q   0)
        (0      1    0      0)
        (sin q  0  cos q    0)
        (0      0    0     1)
Note:
To construct a unique rotation matrix, you will need 4 vectors (not on the same plane) in 3D.


