Let's clarify a number of concepts here to help elucidate how to solve this problem. We know for starters that $$\left[\vec{L}^2,L_i\right]=0,\quad i\in(x,y,z),$$ meaning that we can define simultaneous eigenstates of $\vec{L}^2$ and one of the angular momentum projections $L_i$. Moreover, any other angular momentum projection will also suffice, defined through $$L_{\vec{n}}=\vec{n}\cdot\vec{L}.$$
By spherical symmetry, we do not expect any direction to be privileged, so we expect the set of eigenvalues of $L_{\vec{n}}$ and $\vec{L}^2$ to be independent of the direction $\vec{n}$. That means that we should always be able to define a set of eigenstates $|l,m\rangle_{\vec{n}}$ for any direction such that
$$\vec{L}^2|l,m\rangle_{\vec{n}}=l(l+1)|l,m\rangle_{\vec{n}},\quad L_{\vec{n}}|l,m\rangle_{\vec{n}}=m|l,m\rangle_{\vec{n}},\quad l\in1/2,1,3/2,\cdots,\quad m\in(-l,-l+1,\cdots,l).$$
How do we know this can be done? We define rotation operators $$R(\theta,\vec{u})=\exp(-i\theta L_{\vec{u}}),$$ which [must be unitary due to Wigner's theorem and] have the nice property
$$R\vec{L} R^\dagger=\mathbf{R}^T\vec{L} \qquad\Rightarrow RL_{\vec{n}} R^\dagger=L_{\mathbf{R}\vec{n}},$$ where $\mathbf{R}^T$ is the transpose of a $3\times 3$ rotation matrix that rotates a vector by angle $\theta$ about axis $\vec{u}$. Then, we can do some manipulations starting from the known properties of the states $|l,m\rangle_z$:
\begin{align}
m \left(R|l,m\rangle_z\right)&=R m|l,m\rangle_z\\
&=R L_z|l,m\rangle_z\\
&=R L_zR^\dagger R|l,m\rangle_z\\
&=L_{\mathbf{R}z} \left(R|l,m\rangle_z\right).
\end{align} This tells us that the rotated angular momentum projection operators $L_{\vec{n}}=L_{\mathbf{R}z}$ have the same eigenvalues $m$ before and after the rotation, with the relationships between the eigenstates being given by
$$|l,m\rangle_{\mathbf{R}z}=R|l,m\rangle_z.$$
Given this mathematical machinery, we can always define the eigenstates $|l,m\rangle_{\vec{n}}$ as above, then use the completeness of the eigenstates to write
$$|l,m\rangle_{\vec{n}}=\sum_k |l,k\rangle_{z}\langle l,k|l,m\rangle_{\vec{n}}\equiv \sum_k c_k |l,k\rangle_{z}.$$ In your question, you are looking for these coefficients $c_k=\,_{z}\langle l,k|l,m\rangle_{\vec{n}}=\,_{z}\langle l,k|R|l,m\rangle_{z}$. Does our mathematical machinery help with that? Yes!
These coefficients are exactly the elements of the Wigner $D$-matrices, which have been tabulated and can be calculated analytically in various ways. Your computation helps calculate them for $l=1$.