You’re essentially looking for the difference in the spherical coordinates of the two vectors, so you can use the standard formulas for conversion from Cartesian coordinates: $$r=\sqrt{x^2+y^2+z^2} \\ \theta = \arccos{\frac zr} \\ \phi = \arctan{\frac yx}.$$ Remember that $\theta$ is measured from the positive $z$-axis, and you’ll need to adjust $\Delta\theta$ to be in the range $(-\pi,\pi]$ if you want the use both clockwise and counterclockwise rotations about the $z$-axis for the first rotation.
If all you really need is the two rotations, though, they can be constructed without explicitly computing either angle. The basic idea is to build sets of orthonormal bases that represent the “before” and “after” states. For the first rotation, project $\mathbf u$ and $\mathbf v$ onto the $x$-$y$ plane and normalize to get $\mathbf u' = (u_x^2+u_y^2)^{-1/2}(u_x,u_y,0)^T$ and $\mathbf v' = (v_x^2+v_y^2)^{-1/2}(v_x,v_y,0)^T$, respectively. The required rotation matrix is then $$R_1 = \begin{bmatrix}\mathbf u'&\mathbf e_z\times\mathbf u'&\mathbf e_z\end{bmatrix} \begin{bmatrix}\mathbf v'&\mathbf e_z\times\mathbf v'&\mathbf e_z\end{bmatrix}^{-1} = \begin{bmatrix}\mathbf u'&\mathbf e_z\times\mathbf u'&\mathbf e_z\end{bmatrix} \begin{bmatrix}\mathbf v'^T \\ (\mathbf e_z\times\mathbf v')^T \\ \mathbf e_z^T\end{bmatrix},$$ with $\mathbf e_z = (0,0,1)^T$, the unit $z$-vector. (Note that computing the cross product with $\mathbf e_z$ are a matter of swapping coordinates and negating one of them.)
The second rotation is constructed in a similar manner, except that you’re now rotating $R_1\mathbf v$ onto $\mathbf u$, so you’ll need to normalize both of those vectors and use $\mathbf w = R_1\mathbf v\times\mathbf u$, normalized as the rotation axis.
If needed, the rotation angle $\theta$ can be extracted from each of the resulting rotation matrices using the identity $\operatorname{tr}R = 1+2\cos\theta$.