From the expression of the Hamiltonian you can see that $\hat{H}=\hat{H}_{1}+\hat{H}_{2}$. Also you have the following commutation relation $[\hat{H}_{1},\hat{H}_{2}]=0$. Thus, you can label the eigenstates as $|n_{1}n_{2},SM_{s}\rangle$. Where $|SM_{s}\rangle$ is the two particle total spin state.
Now, the ket corresponding to the ground state is $|n_{1}n_{2},SM_{s}\rangle=|00,00\rangle$. To see where the spin comes into play, we look at the first excited state. The first excited energy is $E_{1}=E(10)=E(01)$ with four possible kets:
For the singlet state $\left(\frac{1}{\sqrt{2}}|10\rangle+\frac{1}{\sqrt{2}}|01\rangle\right)|00\rangle$
For the triplet state $\left(\frac{1}{\sqrt{2}}|10\rangle+\frac{1}{\sqrt{2}}|01\rangle\right)|S=1,M_{s}=0,\pm1\rangle$
Having fermions, the antisymmetric wave function is
$$\psi=\frac{1}{\sqrt{2}}(\psi_{1}(x_{1})\psi_{2}(x_{2})-\psi_{1}(x_{2})\psi_{2}(x_{1}))$$
(there's a plus in your wave function and that is for integer spin particles). This wave function can be split into a spatial and spin part. Being an antisymmetric wave function, when the spatial part is symmetric the spin part is antisymmetric and vice versa.
$$\psi(x_{1},x_{2},M_{1},M_{2})=\left\{
\begin{array}{ll}
\psi^{S}(x_{1},x_{2})\chi^{A}(M_{1},M_{2})\\
\psi^{A}(x_{1},x_{2})\chi^{S}(M_{1},M_{2})
\end{array}\right.$$
where the spatial part is
$$\psi^{S}(x_{1},x_{2})=\frac{1}{\sqrt{2}}(\psi_{1}(x_{1})\psi_{2}(x_{2})+\psi_{1}(x_{2})\psi_{2}(x_{1}))$$
$$\psi^{A}(x_{1},x_{2})=\frac{1}{\sqrt{2}}(\psi_{1}(x_{1})\psi_{2}(x_{2})-\psi_{1}(x_{2})\psi_{2}(x_{1}))$$
and the spin part
$$\chi^{A}(M_{1},M_{2})=\frac{1}{\sqrt{2}}(\uparrow_{1}\downarrow_{2}-\uparrow_{2}\downarrow_{2}) \hspace{5mm} M_{s=}0,\hspace{5mm} S=0$$
$$\chi^{S}(M_{1},M_{2})=\left\{\begin{array}{ll}
\uparrow_{1}\uparrow_{2},\hspace{29mm} M_{s}=1\\
\frac{1}{\sqrt{2}}(\uparrow_{1}\downarrow_{2}+\uparrow_{2}\downarrow_{2}) \hspace{5mm} M_{s=0}\\
\downarrow_{1}\downarrow_{2} \hspace{29mm} M_{s}=-1
\end{array}\right \} , S=1$$
It is true that you obtain an equation for $u$ that is exactly equal to the equation for the 1D harmonic oscillator; what changes are boundary conditions.
In fact, in solving these kind of equations, you require that the radial solution $R(r)$ goes like a certain power, that turns out to be
$$
R(r) \xrightarrow[r \to 0]{} r^l\, .
$$
Then necessarily we have $u(r)\xrightarrow[r \to 0]{} r^{l+1}$, that in particular for $l=0$ means
$$
u(0)=0\,, \qquad u(r)\xrightarrow[r \to 0]{} r \, \,.
$$
From here it can be noticed that this solution does not correspond to the ground state of the 1D harmonic oscillator, that being a Gaussian is not null at $r=0$.
The first eigenfunction of the 1D harmonic oscillator that fulfills the boundary conditions is actually the one corresponding to $n=1$, with an energy given by
$$
E=\hbar \omega\left(1+\frac12\right)= \frac32 \hbar \omega \, ,
$$
which is the expected result.
Best Answer
You can definitely represent a 3d QHO wavefunction as a composition of radial components and angular components (spherical harmonics).