You wrote down a matrix that maps $(0, 0, 1)$ to $\vec x$. That is the inverse of the matrix that maps $\vec x$ to $(0, 0, 1)$.
Edit: it's possible I misunderstood, so let's go through it.
You want to define a set of cartesian axes $e_i$ such that $\vec x = e_z$. You know $\vec x = {x'}^i {e'}_i$ in an old basis that you're rotating from. The bases are related by a rotation operator $\underline R({e'_i}) = e_i$.
So, let's look at the components of ${e'}_z$ in the new basis. We want to evaluate ${e'}_z \cdot e_j = {e'}_z \cdot \underline R({e'}_j) = \underline R^{-1}({e'}_z) \cdot {e'}_j$. This means that to find the components of ${e'}_z$ on the new basis, you can equivalently inverse rotate it and find the components on the old basis. But $\underline R^{-1}({e'}_z) \neq \vec x$--picture this if you don't follow. That you got the erroneous result that they were the same tells me your matrix is not correct. It is not what you thought it was; it's probably assuming you want to rotate vectors and not the axes or some other similar, subtle disconnect.
Edit edit: actually, I'm not sure I understand what you're saying you found ${e_z}'$ to be in the new basis. The notation does not make sense.
3rd edit: let's walk through a whole cloth calculation of the rotation operator. Doing this with two rotations chained is somewhat onerous, if straightforward. I'll use quaternions (spinors, rotors) in a geometric (clifford) algebra.
The idea here is simple: any rotation can be broken down into a double-sided action of a spinor.
$$\underline R(a) = \psi a \psi^{-1}$$
where the products here are clifford products. The spinor $\psi$ is formed as an exponential of a bivector. Given two orthonormal vectors $u, v$ that span the rotation plane, the spinor that rotates through an angle $\theta$ is $\psi = \exp(-u \wedge v \, \theta/2)$.
First, we find the vector that lies in the plane $\vec x \wedge e_z$ that is orthonormal to $e_z$. The expression is rather simple:
$$\rho = (\vec x \wedge e_z) \cdot e^z = (\sin \theta \cos \phi e_{xz} + \sin \theta \sin \phi e_{yz}) \cdot e^z = (\sin \theta \cos \phi e_x + \sin \theta \sin \phi e_y)$$
Exactly as you'd expect. This vector has magnitude $\sin \theta$, so the normalized vector is
$$\hat \rho = e_x \cos \phi + e_y \sin \phi$$
Because we've decomposed the problem into spherical coordinates, we know the angle between $\vec x$ and $e_z$ is $\theta$. We can then write the spinor as
$$\psi = \exp(-\hat \rho e_z \theta/2) = \cos \frac{\theta}{2} - \hat \rho e_z \sin \frac{\theta}{2}$$
With $\hat \rho e_z = e_{xz} \cos \phi + e_{yz} \sin \phi$. Now let's get down to evaluating the operator's components.
$$\begin{align*}\underline R(e_x) &= \Big ( \cos \frac{\theta}{2} - \hat \rho e_z \sin \frac{\theta}{2} \Big) e_x \Big ( \cos \frac{\theta}{2} + \hat \rho e_z \sin \frac{\theta}{2} \Big) \\
&= e_x \cos^2 \frac{\theta}{2} + e_z \sin \theta \cos \phi - \sin^2 \frac{\theta}{2} (e_x \cos 2\phi + e_y \sin 2\phi) \end{align*}$$
The other basis vectors are mapped to
$$\begin{align*} \underline R(e_y) &= e_y \cos^2 \frac{\theta}{2} + e_z \sin \theta \sin \phi + \sin^2 \frac{\theta}{2} (e_y \cos 2\phi - e_x \sin 2\phi) \\ \underline R(e_z) &= e_z \cos \theta - (e_x \cos \phi + e_y \sin \phi) \sin \theta \end{align*}$$
Checking that this is correct for the full vector is a bit involved, but I've checked a few important limits ($\phi = 0, \pi/2$) and the behavior seems to check out.
