Let $v$ be your vector. If $x=y=0$, then, since it is an unit vector, $z=\pm1$. If $z=1$, you can take the identity, and if $z=-1$, you can take a half-turn around any vector orthogonal to $(0,0,1)$.
Now, I will assume that $(x,y)\ne(0,0)$. Let$$w_1=\left(\frac y{\sqrt{x^2+y^2}},\frac{-x}{\sqrt{x^2+y^2}},0\right),$$which is an unit vector which is orthogonal to both $v$ and $(0,0,1)$, and let$$w_2=\left(\frac x{\sqrt{x^2+y^2}},\frac y{\sqrt{x^2+y^2}},0\right),$$which is also an unit vector and which is orthogonal to both $w_1$ and $(0,0,1)$. So, $\bigl\{w_1,w_2,(0,0,1)\bigr\}$ is an orthonormal basis of $\Bbb R^3$. Take $\theta\in[0,\pi]$ such that $\cos(\theta)=z$ and that $\sin(\theta)\geqslant0$ (which implies that $\sin(\theta)=\sqrt{1-z^2}=\sqrt{x^2+y^2}$) and note that, since $\cos\theta=z=\bigl\langle(v,(0,0,1)\bigr\rangle$, $\theta$ is the angle between $v$ and $(0,0,1)$. Consider the rotation $R$ around $w_1$ with angle $\theta$. More precisely, consider the linear map $R\colon\Bbb R^3\longrightarrow\Bbb R^3$ such that:
- $R(w_1)=w_1$;
- $R(w_2)=\cos(\theta)w_2-\sin(\theta)(0,0,1)$;
- $R(0,0,1)=\sin(\theta)w_2+\cos(\theta)(0,0,1)=v$.
Then $R$ is a rotation that maps $(0,0,1)$ into $v$. The matrix of $R$ with respect to the standard basis of $\Bbb R^3$ is$$\begin{bmatrix}z+\frac{y^2}{z+1} & -\frac{x y}{z+1} & x \\ -\frac{x y}{z+1} & z+\frac{x^2}{z+1} & y \\ -x & -y & z\end{bmatrix}.$$What you're after is $R^{-1}$, which is the transpose of $R$, since $R$ is an orthogonal matrix. So, an answer to your question is$$\begin{bmatrix}z+\frac{y^2}{z+1}&-\frac{x y}{z+1}&-x\\-\frac{x y}{z+1}&z+\frac{x^2}{z+1}&-y\\x&y&z\end{bmatrix}.$$