We introduce the ladder operators $a^\dagger_i, a_i$ such that \begin{align}x_i & = \sqrt{\frac{\hbar}{2m\omega}} (a^\dagger_i + a_i) \\
p_i & = i\sqrt{\frac{\hbar m\omega}{2}} (a^\dagger_i - a_i) \end{align}
where $i = 1,2,3$. The commutators are of course $$[a_i, a^\dagger_j] = \delta_{ij}.$$
Then the angular momentum operator is $$L_i = \epsilon_{ijk}x_j p_k$$
with $\epsilon_{ijk}$ the Levi-Civita symbol and sums over $j, k$ implied. On expanding $x_j p_k$ only $-a^\dagger_j a_k$ and $a_j a^\dagger_k$ contribute, since $a_k a_j$ is symmetric in $k,j$. These two terms give equal contributions since their commutator is symmetric in $k,j$. It follows that $$L_i = -i\hbar\epsilon_{ijk} a^\dagger_j a_k.$$
Now defining $$a_+ = \frac{-1}{\sqrt 2} (a_x - ia_y) \qquad a_- = \frac{1}{\sqrt 2}(a_x + ia_y)$$
we have $[a_\pm, a_\pm^\dagger]= 1$ and $$L_z = \hbar(a^\dagger_+ a_+ - a^\dagger_- a_-).$$
It is rather clear that $a^\dagger_\pm$ raise $N$ by $1$, and $a_\pm$ adds an excitation with $L_z = \pm\hbar$: $L_z$ is the difference between the number operators corresponding to $a_\pm$.
Using these operators you can in principle work out the matrix for $L_z$ (and also $L_x$ and $L_y$) and $L^2$. Since the $L_i$ operators contain only products one creation and one annihilation operator they do not connect states with different $N$. It follows that neither does $L^2$, so you can consider each $N$ separately. Once you have these matrices, diagonalizing them will tell you how to express the $|N, l, m\rangle$ in terms of the $|n_x,n_y, n_z\rangle$.
Note that for each $N$ there are $(N+2)(N+1)/2$ states, so you probably don't want to do this by hand except maybe for $N = 2$ (the $N = 0,1$ cases are trivial). Maybe you can get Mathematica or Maple to do it for you for some larger $N$.
For $N = 1$ we can calculate as follows. $$a_+ |0,0,1\rangle =\frac{-1}{\sqrt 2}(a_x - ia_y)|0,0,1\rangle = 0.$$
What we have used here is that $$a_x |n_x, n_y, n_z\rangle = \sqrt{n_x}|n_x-1,n_y,n_z\rangle$$ and similarly for $y,z$. We obtain the same with $a_-$, so this means that $L_z |0,0,1\rangle = 0$, so the state $|0,0,1\rangle$ has $m = 0$. For $|1,0,0\rangle$, we have $$a_+|1,0,0\rangle = -a_- |0,1,0\rangle = -\frac{1}{\sqrt 2} |0,0,0\rangle$$
from which we get $$a^\dagger_+ a_+ |1,0,0\rangle = \frac{1}{2} \big( |1,0,0\rangle + i |0,1,0\rangle\big)$$ $$a^\dagger_- a_- |1,0,0\rangle = \frac{1}{2} \big(|1,0,0\rangle - i|0,1,0\rangle\big).$$ Thus $$L_z |1,0,0\rangle = i\hbar |0,1,0\rangle.$$ Now you can probably work out on your own that $L_z |0,1,0\rangle = -i\hbar |1,0,0\rangle.$ This gives the matrix for $L_z$ on $N = 1$ states as $$L_z = \begin{pmatrix}0 & i & 0 \\ -i & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}.$$ The eigenvalues of $L_z$ are given by the solutions to $$m(m^2-1) = 0$$ which are $m = 0, \pm 1$. We already know that $|0,0,1\rangle$ is the eigenvector corresponding to $m = 0$. To find the eigenvector corresponding to $m = \pm 1$, we have to solve the system of linear equations \begin{align}iy & = \pm x \\ -ix & = \pm y\end{align} which just says $x = \pm i y$ (the two equations are equivalent). Thus $$L_z \big( |1,0,0\rangle \pm i |0,1,0\rangle\big) = \pm\hbar \big( |1,0,0\rangle \pm i |0,1,0\rangle\big).$$
Since we found three states, with $m = -1,0,1$, we must have $l = 1$.
Of course in this simple case we could also have reasoned like this: $a^\dagger_\pm$ adds an excitation with angular momentum $\pm \hbar$. Knowing that $L_z|0,0,0\rangle = 0$, we get states with $m = \pm 1$ just by acting with $a^\dagger_\pm$ on $|0,0,0\rangle$. Indeed up to a normalization this is just what we found.
