As an alternative to your method, which works, you can use spherical coordinates to express the axis of rotation $\mathbf{v}$ as follows
$ \mathbf{v} = [ \sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta ]^T $
Then construct the following two unit vectors that are orthogonal to $\mathbf{v}$ and also to each other,
$ \mathbf{u_1} = [ \cos \theta \cos \phi, \cos \theta \sin \phi, -\sin \theta ]^T$
$ \mathbf{u_2} = [-\sin \phi, \cos \phi, 0 ]^T $
Now assemble the three vectors together into one matrix $R$ as follows
$ R = [ \mathbf{v} , \mathbf{u_1}, \mathbf{u_2} ] $
The columns of matrix $R$ represent a basis for $\mathbb{R}^3$. In this basis, the coordinate of a given point $P = (x, y, z)$ is
$Q =(x',y',z') = R^T \ P $
Now we want to rotate $P$ about $\mathbf{v}$ by an angle $\lambda$. This means that we're rotating $Q$ about the the $x'$ axis (the first column of $R$). The rotated point will have its coordinates (in the basis $R$) as
$ Q' = R_{x'} \ Q = \begin{bmatrix} 1 && 0 && 0 \\ 0 && \cos \lambda && - \sin \lambda \\ 0 && \sin \lambda && \cos \lambda \end{bmatrix} \begin{bmatrix} x' \\ y' \\ z' \end{bmatrix} $
Finally express the rotated point in the standard basis, as follows
$ P' = R \ Q' = R \ R_{x'} \ Q = R \ R_{x'} \ R^T \ P $
So that the overall rotation matrix in the standard basis is given by
$ R_{v} = R \ R_{x'} \ R^T $
This expression requires the evaluation of the matrix $R$ with $\mathbf{v}, \mathbf{u_1} , \mathbf{u_2} $ as its columns. In what follows, I'll eliminate the requirement to calculate $\mathbf{u_1}$ and $\mathbf{u_2}$, and I will produce an expression that depends only on $\mathbf{v}$ and the angle $\lambda$.
Recall that $R = [ {\mathbf{v} , \mathbf{u_1}, \mathbf{u_2}} ] $, so the above expression for the overall rotation matrix becomes,
$R_{v} = {\mathbf{v} \mathbf{v}}^T + \cos(\lambda) \bigg( {\mathbf{u_1} \mathbf{u_1}}^T + {\mathbf{u_2} \mathbf{u_2}}^T \bigg) + \sin(\lambda) \bigg( - {\mathbf{u_1} \mathbf{u_2}}^T + {\mathbf{u_2} \mathbf{u_1}}^T \bigg) $
Since ${R R}^T = I$, then ${\mathbf{u_1} \mathbf{u_1}}^T + {\mathbf{u_2} \mathbf{u_2}}^T = I - {\mathbf{v} \mathbf{v}}^T $
The sine term needs some careful consideration. If $P$ is any point in $\mathbb{R}^3 $, then its coordinates in the basis defined by the columns of $R$ is
$ Q = (x',y', z') = R^T P $
i.e. $ P = x' \mathbf{v} + y' \mathbf{u_1} + z' \mathbf{u_2} $
where
$ x' = \mathbf{v}^T P $
$ y' = \mathbf{u_1}^T P $
$ z' = \mathbf{u_2}^T P $
The vector $W = \bigg( - {\mathbf{u_1} \mathbf{u_2}}^T + {\mathbf{u_2} \mathbf{u_1}}^T \bigg) P = y' \mathbf{u_2} - z' \mathbf{u_1} $
Now consider the cross product:
$\mathbf{v} \times P = \mathbf{v} \times (x' \mathbf{v} + y' \mathbf{u_1} + \mathbf{u_2} ) = 0 + y' \mathbf{u_2} - z' \mathbf{u_1} $
Comparing the two expressions above, we deduce that
$\bigg( - {\mathbf{u_1} \mathbf{u_2}}^T + {\mathbf{u_2} \mathbf{u_1}}^T \bigg) P = \mathbf{v} \times P $
This cross product between $\mathbf{v}$ and $P$ can be expressed as a matrix multiplication as follows
$ \mathbf{v} \times P = S_v P = \begin{bmatrix} 0 && - v_z && v_y \\ v_z && 0 && - v_x \\ - v_y && v_x && 0 \end{bmatrix} P $
In summary, the rotation matrix is given by
$ R_v(\lambda) = \mathbf{v} \mathbf{v}^T + \cos(\lambda) (I - \mathbf{v} \mathbf{v}^T)+ \sin(\lambda) S_v $
And this is the well-known formula for the rotation matrix known as Rodrigues' Rotation Matrix Formula.
Best Answer
For the Z axis rotation, you want to rotate through an angle whose tangent is $\frac{x}{y}$, so the rotation matrix is $$ \langle x,y,z\rangle \cdot \begin{pmatrix} \frac{y}{\sqrt{x^2+y^2}}& \frac{x}{\sqrt{x^2+y^2}}& 0 \\ -\frac{x}{\sqrt{x^2+y^2}} & \frac{y}{\sqrt{x^2+y^2}} & 0 \\ 0 & 0 & 1 \end{pmatrix} = \langle 0, \sqrt{x^2+y^2}, z\rangle.$$ For the second rotation, the rotation around the $x$-axis is clockwise through the complement of the angle whose tangent is $\frac{\sqrt{x^2+y^2}}{z}$. For a general point $\langle 0, y, z\rangle$, the associated rotation matrix is $$\begin{pmatrix} 1 & 0 & 0 \\ 0 & -\frac{z}{\sqrt{y^2+z^2}} & -\frac{y}{\sqrt{y^2+z^2}} \\ 0 & \frac{y}{\sqrt{y^2+z^2}} & \frac{z}{\sqrt{y^2+z^2}} \end{pmatrix}.$$ Substituting $\sqrt{x^2+y^2}$ for $y$ in this matrix gives for the $x$-axis rotation matrix $$\begin{pmatrix} 1 & 0 & 0 \\ 0 & -\frac{z}{\sqrt{x^2+y^2+z^2}} & -\frac{\sqrt{x^2+y^2}}{\sqrt{x^2+y^2+z^2}} \\ 0 & \frac{\sqrt{x^2+y^2}}{\sqrt{x^2+y^2+z^2}} & \frac{z}{\sqrt{x^2+y^2+z^2}} \end{pmatrix}.$$ Multiplying $\langle x, y, z\rangle$ by these matrices in turn gives the vector $\langle 0, 0, \frac{-x^2 - y^2 + z^2}{\sqrt{x^2 + y^2 + z^2}}\rangle$.