Absorbing the irrelevant ħω constants into the normalization of the suitable quantities, for the 3D isotropic oscillator, $\epsilon=n+3/2$, while for each n the degeneracy is $(n+1)(n+2)/2$; (see SE ). Scoping the power behavior of a large quasi-continuous n, leads you to the answer.
The number of states then goes like $N\propto n^3 \propto \epsilon^3$, and hence the density of states like $dN/d\epsilon\propto \epsilon^2$.
Let us start from the one-dimensional case. The the parity operator $R$
is defined as $$(R\psi)(x) := \psi(-x)$$ for every wavefunction $\psi \in L^2(\mathbb R, dx)$.
$R$ is selfadjoint and unitary and, directly from its definition, we get $RX=-XR$ and $RP=-PR$ where $X$ and $P$ are the position and momentum operators.
Since the creation operator $a^\dagger$ and the annihilation operator $a$ of the harmonic oscillator are linear combinations of $X$ and $P$, the same commutation relations are valid for $a$ and $a^\dagger$.
In particular, $$Ra^\dagger=-a^\dagger R\:.\tag{1}$$
Since $R|0\rangle =|0\rangle$ as it arises by direct inspection in position representation (the wavefunction of the ground state is proportional to $e^{-cx^2}$ which is an even function), we have from (1) that $$R|n\rangle = R \frac{(a^\dagger)^n}{\sqrt{n!}}|0\rangle = (-1)^n \frac{(a^\dagger)^n}{\sqrt{n!}}R|0\rangle = (-1)^n |n\rangle\:.\tag{2}$$
In the considered case $$|n_1 n_2 n_3\rangle = |n_1\rangle\otimes |n_2\rangle \otimes |n_3\rangle\tag{3}$$ and the parity operator $R$ is the tensor product of the corresponding three parity operators $$R = R_1 \otimes R_2 \otimes R_3$$ so that the parity of the state $|n_1 n_2 n_3\rangle$ is $(-1)^{n_1+ n_2 +n_3}$. Indeed, from (1),(2), and (3):
$$R |n_1 n_2 n_3\rangle = (-1)^{n_1+ n_2 +n_3} |n_1 n_2 n_3\rangle\:.$$
This is the parity of the eigenvalue, unless there is accidental degeneracy: $$E_{n_1, n_2, n_3} = E_{m_1, m_2, m_3}\quad \mbox{for some $(n_1, n_2, n_3) \neq (m_1, m_2, m_3)$,}$$ in that case the eigenvalue may have not defined parity. Accidental degeneracy arises under swap of $n_1$ and $n_2$ in the considered case, however the parity is invariant under this operation.
Best Answer
What you have here could be described as a subset sum problem. Given $n$ can take any integer value (not including zero), you have the set of squares up to $36$,
$S = \{1,4,9,16,25,36\}$
and you wish to find subsets of three which sum to $41$. Looking at the subset sum problem this can not be solved analytically but algorithms can be employed.
To do this vigorously you need to form a 3D matrix $6 \times 6 \times 6$ by summing together the relevant squares and then read off the indices where $41$ is achieved, a task made easier with a bit easier with a computer script.
You must also note that any summation where $n_x \neq n_y \neq n_z$ will form a 6-fold permutation, sets with 2 distinct values will form a 2-fold permutation, whereas for $n_x = n_y = n_z$ the solution is unique.