One approach is to use quaternions.
Construct the unit vector along the rotation axis. Construct a unit quaternion using that vector and half the desired rotation. If the unit vector is u, and theta is half the rotation angle, the quaternion will be cos(theta)+sin(theta)*u. Call it Q. Then you can rotate andy vector V by extending it to a quaternion (with a zero scalar part) and calculating Q*V*Q^-1. Note that for unit quaternions, the reciprocal (Q^-1) is the same as the conjugate.
With the normal coordinate system, this will be a counterclockwise rotation as viewed looking down the direction vector. If you want to consider a rotation of the axes and get the representation of the same vectors in the new coordinate system, you just take the negative of the angle. You can also work this out symbolically for each of the three basis vectors in order to calculate the rotation matrix.
Tom Gutman