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.


