How to calculate the angular momentum states of isotropic quantum harmonic oscillator?
As you have probably already worked out, the two eigenspaces that you're interested in, $N=2$ and $N=3$, are formed by a direct sum of subspaces of different angular momentum; thus, the $N=2$ eigenspace has one $l=2$ and one $l=0$ subspace, of dimensions $5+1=6$, and the $N=3$ eigenspace has one $l=3$ and one $l=1$ subspace, of dimensions $7+3=10$. You therefore have two distinct tasks: moving around within each subspace, and jumping to a different representation.
The first task is relatively easy, and in fact you do not need much knowledge of the internal structure of the hamiltonian to do this. You already know that $H$ commutes with $\mathbf L$, which guarantees you the shared eigenbasis $|nlm⟩$, but more than that you know that acting on the $|nlm⟩$ with components in the angular momentum, and particularly the ladder operators $L_\pm = L_x+iL_y$, will keep you inside that subspace. In this vein, for example, you can take the $|3,3,3⟩$ state you've already found to get $$L_-|3,3,3⟩\propto |3,3,2⟩,$$ and so on.
What you really need to implement this in your language, though, is to bring the $L_\pm$ into the language of your creation and annihilation operators. We've already done most of the heavy lifting in the previous question, which gave us the identity $$ L_i=-i\hbar\varepsilon_{ijk}a_j^\dagger a_k^\phantom{\dagger}, $$ and now we just need to apply this to the ladder operators: \begin{align} L_\pm & = L_x \pm i L_y \\ & = -i\hbar\varepsilon_{1jk}a_j^\dagger a_k^\phantom{\dagger} \pm \hbar\varepsilon_{2jk}a_j^\dagger a_k^\phantom{\dagger} \\ & = -i\hbar(a_2^\dagger a_3^\phantom{\dagger} -a_3^\dagger a_2^\phantom{\dagger} )\pm \hbar(a_3^\dagger a_1^\phantom{\dagger}-a_1^\dagger a_3^\phantom{\dagger}) \\ & = \hbar(\mp a_1^\dagger-ia_2^\dagger)a_3^\phantom{\dagger} + \hbar a_3^\dagger(\pm a_1^\phantom{\dagger}+ia_2^\phantom{\dagger}) \\ & = \sqrt{2}\hbar(a_\pm^\dagger a_0^\phantom{\dagger} + a_0^\dagger a_\mp^\phantom{\dagger}), \end{align} where I've used the forms $a_\pm = \frac{1}{\sqrt{2}}(\mp a_1 +ia_2)$ and $a_\pm^\dagger =\frac{1}{\sqrt{2}}(\mp a_1^\dagger -ia_2^\dagger) $. This result for $L_\pm$ can be understood rather easily: you destroy one quantum with $a_0$ and then recreate it with $\pm$more angular momentum, or you destroy one quantum with $a_\mp$ and then recreate it with $\pm1$ more unit of angular momentum. Easy!
As applied to your $N=3$ state, $|3,3,3⟩ = \frac{1}{\sqrt{3!}}(a_+^\dagger)^3|0⟩$, you can then proceed to get \begin{align} |3,3,2⟩ &= \frac{1}{\sqrt{6}}L_- |3,3,3⟩ \\ & = \frac{1}{\sqrt{6}} \, \sqrt{2} \hbar(a_-^\dagger a_0^\phantom{\dagger} + a_0^\dagger a_+^\phantom{\dagger}) \ \frac{1}{\sqrt{3!}}(a_+^\dagger)^3|0⟩ \\ & = \frac{1}{3\sqrt{2}}\hbar(0 + a_0^\dagger a_+^\phantom{\dagger}) \cdot (a_+^\dagger)^3|0⟩ \\ & = \frac{1}{3\sqrt{2}}\hbar\, a_0^\dagger \left((a_+^\dagger)^3 a_+^\phantom{\dagger}+3(a_+^\dagger)^2\right) |0⟩ \\ & = \frac{1}{\sqrt{2!}} \hbar \, a_0^\dagger (a_+^\dagger)^2 |0⟩, \end{align} and so on, until you reach $|3,3,-3⟩ = \frac{1}{\sqrt{3!}}(a_-^\dagger)^3|0⟩$.
The other subspace, with $l<N$, is a bit more tricky, because you cannot get it from the naive argument of just doing $(a_+^\dagger)^{N} |0⟩$ and then using the angular momentum algebra to cover you.
I always find the emergence of these shells to be easier to understand by looking at the way the cartesian components combine, so let me start by building the $N=2$, $l=0$ state directly on the cartesian basis $|n_x,n_y,n_z⟩_\mathrm{C}$. In particular, consider the state $$ |\psi⟩ = \frac{1}{\sqrt{3}} \left( |2,0,0⟩_\mathrm{C} + |0,2,0⟩_\mathrm{C} + |0,0,2⟩_\mathrm{C}\right), $$ which has a position-space representation \begin{align} ⟨x,y,z|\psi⟩ & = \frac{1}{\sqrt{3}} \left( ⟨x,y,z|2,0,0⟩_\mathrm{C} + ⟨x,y,z|0,2,0⟩_\mathrm{C} + ⟨x,y,z|0,0,2⟩_\mathrm{C}\right) \\ & = \frac{1}{\sqrt{3}} \bigg[ H_2(x)+H_2(y)+H_2(z)\bigg]e^{-r^2/2} \\ & \propto \frac{1}{\sqrt{3}} \bigg[ (4x^2-2)+(4y^2-2)+(4z^2-2)\bigg]e^{-r^2/2} \\ & = \frac{1}{\sqrt{3}} \left(4r^2-6\right)e^{-r^2/2} , \end{align} so it depends only on $r$ and therefore belongs to the $l=0$ representation. This is all well and good, but how do we represent this using the bosonic operators? This can be read directly from the state above, by rephrasing it as \begin{align} |2,0,0⟩ & = \frac{1}{\sqrt{3!}}\left((a_1^\dagger)^2+(a_2^\dagger)^2+(a_3^\dagger)^2\right)|0⟩ \\ & = \frac{1}{\sqrt{3!}}\left(\frac12(-a_+^\dagger +a_-^\dagger)^2-\frac12(a_+^\dagger +a_-^\dagger)^2+(a_0^\dagger)^2\right)|0⟩ \\ & = \frac{1}{\sqrt{3!}}\left((a_0^\dagger)^2-2a_+^\dagger a_-^\dagger\right)|0⟩ \end{align} where I've used $a_1^\dagger = \frac{1}{\sqrt{2}}\left(-a_+^\dagger +a_-^\dagger\right)$ and $a_2^\dagger = \frac{i}{\sqrt{2}}\left(a_+^\dagger +a_-^\dagger\right)$. This makes a lot of sense: $|2,0,0⟩$ is generated up from the vacuum state through two pathways that increase $N$ by two but leave $L_0$ unchanged, $(a_0^\dagger)^2$ and $a_+^\dagger a_-^\dagger$, and then those two are superposed in a way that will make $L^2$ vanish. There is an obvious question, though - what is that operator, exactly, and how do we generalize it?
The answer here is that the above example is illustrative, but it is not quite the right approach, yet. We don't really want operators that will take us from one energy eigenspace to another, because those will not commute with the hamiltonian and will therefore have some relatively complicated algebraic properties. Instead, what we really want is a way to generate the $|2,0,0⟩$ state by jumping down the $l$ ladder from the $N=2$, $l=2$ space: if we do this, then our new class of ladder operators can commute with the hamiltonian (but not with $\mathbf L$).
This gives you a couple of clear candidates, because the isotropic oscillator has very few independent symmetries: the angular momentum, which we've already used, the Laplace-Runge-Lenz vector, and something called the Fradkin tensor, which is defined as $$ F_{ij} = p_ip_j +x_ix_j $$ in the cartesian frame, and which turns out to commute with the hamiltonian (which you should check explicitly). The Fradkin tensor tends to play nicer with the harmonic oscillator than the Runge-Lenz vector, and I've tried to make it mesh below but to be honest it's not quite there so you're going to need to fill in some blanks.
To bring this over to our old language, let's start by changing the quadratures over into bosonic operators, giving us \begin{align} F_{ij} & = p_ip_j +x_ix_j \\ & = \frac{a_i^\phantom{\dagger}-a_i^\dagger}{\sqrt{2}\,i} \frac{a_j^\phantom{\dagger}-a_j^\dagger}{\sqrt{2}\, i} + \frac{a_i^\phantom{\dagger}+a_i^\dagger}{\sqrt{2}} \frac{a_j^\phantom{\dagger}+a_j^\dagger}{\sqrt{2}} \\ & = \frac12 \left[-(a_i^\phantom{\dagger}-a_i^\dagger)(a_j^\phantom{\dagger}-a_j^\dagger)+(a_i^\phantom{\dagger}+a_i^\dagger)(a_j^\phantom{\dagger}+a_j^\dagger)\right] \\ & = a_i^\phantom{\dagger}a_j^\dagger+a_i^\dagger a_j^\phantom{\dagger} . \end{align} (Exercise: is this explicitly hermitian? If not, why is it hermitian?) This is, of course, on a cartesian basis for the tensor part, but of we want to be talking about spherical tensors here, so we really want to be using its spherical components, $F_{0,0}$ and $F_{2,m}$; here the scalar component is easy, since $$ F_{0,0} = \sum_{i=1}^3 F_{ii} = 2H $$ is just the hamiltonian, and to get the quadrupole components the easiest route is by analogy with the electric quadrupole $Q_{ij}=x_ix_j$, so you get \begin{align} F_{2,2} &= F_{11}+2iF_{12}-F_{22} \\ &= (a_1^\phantom{\dagger}a_1^\dagger+a_1^\dagger a_1^\phantom{\dagger}) +2i (a_1^\phantom{\dagger}a_2^\dagger+a_1^\dagger a_2^\phantom{\dagger}) - (a_2^\phantom{\dagger}a_2^\dagger+a_2^\dagger a_2^\phantom{\dagger}) \\ & = (a_1^\phantom{\dagger} + i a_2^\phantom{\dagger})(a_1^\dagger+ia_2^\dagger) +(a_1^\dagger+ia_2^\dagger)(a_1^\phantom{\dagger} + i a_2^\phantom{\dagger}) \\ & = -4\,a_+^\dagger a_-^\phantom{\dagger}, \\[2mm] F_{2,1} &= F_{13}+iF_{23} \\ & = (a_1^\phantom{\dagger}a_3^\dagger+a_1^\dagger a_3^\phantom{\dagger}) +i(a_2^\phantom{\dagger}a_3^\dagger+a_2^\dagger a_3^\phantom{\dagger}) \\ & = (a_1^\phantom{\dagger}+ia_2^\phantom{\dagger})a_3^\dagger+(a_1^\dagger +ia_2^\dagger )a_3^\phantom{\dagger} \\ & = \sqrt{2}\left(a_-^\phantom{\dagger}a_0^\dagger-a_+^\dagger a_0^\phantom{\dagger} \right), \\[2mm] F_{2,0} & = 2F_{33}-F_{11}-F_{22} \\ & = 2(a_3^\phantom{\dagger}a_3^\dagger+a_3^\dagger a_3^\phantom{\dagger}) -(a_1^\phantom{\dagger}a_1^\dagger+a_1^\dagger a_1^\phantom{\dagger}) -(a_2^\phantom{\dagger}a_2^\dagger+a_2^\dagger a_2^\phantom{\dagger}) \\ & = 2(a_0^\phantom{\dagger}a_0^\dagger+a_0^\dagger a_0^\phantom{\dagger}) -(a_+^\phantom{\dagger}a_+^\dagger+a_+^\dagger a_+^\phantom{\dagger}) - (a_-^\phantom{\dagger}a_-^\dagger+a_-^\dagger a_-^\phantom{\dagger}) \\ & = 2a_0^\dagger a_0^\phantom{\dagger} -a_+^\dagger a_+^\phantom{\dagger} - a_-^\dagger a_-^\phantom{\dagger} , \\[2mm] \text{and similarly} \\[2mm] F_{2,-1} &= F_{13}-iF_{23} \\ & = (a_1^\phantom{\dagger}a_3^\dagger+a_1^\dagger a_3^\phantom{\dagger}) -i(a_2^\phantom{\dagger}a_3^\dagger+a_2^\dagger a_3^\phantom{\dagger}) \\ & = (a_1^\phantom{\dagger}-ia_2^\phantom{\dagger})a_3^\dagger+(a_1^\dagger -ia_2^\dagger )a_3^\phantom{\dagger} \\ & = \sqrt{2}\left(a_-^\dagger a_0^\phantom{\dagger}-a_+^\phantom{\dagger}a_0^\dagger\right) \quad\text{and} \\[2mm] F_{2,-2} &= F_{11}-2iF_{12}-F_{22} \\ &= (a_1^\phantom{\dagger}a_1^\dagger+a_1^\dagger a_1^\phantom{\dagger}) -2i (a_1^\phantom{\dagger}a_2^\dagger+a_1^\dagger a_2^\phantom{\dagger}) - (a_2^\phantom{\dagger}a_2^\dagger+a_2^\dagger a_2^\phantom{\dagger}) \\ & = (a_1^\phantom{\dagger} - i a_2^\phantom{\dagger})(a_1^\dagger-ia_2^\dagger) +(a_1^\dagger-ia_2^\dagger)(a_1^\phantom{\dagger} - i a_2^\phantom{\dagger}) \\ & = -4\,a_-^\dagger a_+^\phantom{\dagger}. \end{align} These operators obviously commute with the hamiltonian, and while not hermitian they obey $F_{2,m}^\dagger = F_{2,-m}$.
This is as far as I'm going to take it for now, as I'm out of time. However, this should give you a good flavourful of where to take this to complete all the relationships in the relevant algebras; the correct ladder operators are there somewhere and it should just be a matter of whittling things down to get the correct relationships.