Eigenfunctions and energies of a quadratic-bottom finite square well
I see the fatal flaw here. Inside the well, your solution must be a solution to the HO, BUT there are two solutions to the HO. There is the one we're all familiar with, a hermite polynomial times a gaussian, but $e^{x^2}$ is also a valid solution. For the HO, we discard that solution because it isn't normalizeable at the boundaries. In the finite well that's no longer true, your wavefunction can be quite large inside the well, it will still be exponentially decaying outside the well. So a linear combination of these two solutions will be the appropriate solution, where you must of course choose how much of each by matching the wavefunction and it's derivative at both ends of the well to an exponential decay with finite energy E.
- Can we write down the solutions to this potential in closed form?
Yes, you can, though they will likely depend on the eigenvalue, and this will likely be the solution of a transcendental equation that can only be solved numerically. (Note that there's nothing strange about that $-$ you get the same behaviour with the finite square well.)
- Why aren't they the same as in the simple harmonic oscillator case?
Because the Hamiltonian is different. There is absolutely no reason to expect the eigenfunctions to remain the same.
On to the details.
As you've pointed out, this problem is best approached via explicit differential equations, which read $$ -\frac12 \frac{\mathrm d^2}{\mathrm dx^2}\psi + \frac12 x^2 \psi = E\psi \tag{1a} $$ for $|x|<a$ and $$ -\frac12 \frac{\mathrm d^2}{\mathrm dx^2}\psi + V_0 \psi = E\psi \tag{1b} $$ for $|x|>a$. Note that I've dropped a bunch of constants, which are completely unnecessary and which only work to make the algebra harder; in effect, I'm measuring $V_0$ in units of $\hbar\omega$ (of your notation) and $a$ in terms of the ground-state width of the parabolic potential.
Moreover, since the potential is parity-symmetric, the bound eigenfunctions must all be either parity even or parity odd. (For continuum states, which I won't consider here, they can always be chosen that way.) This gives us clear boundary conditions, which are either $\psi(0)=0$ or $\psi'(0)=0$, alternating with the energy level, and $\psi(x)\to 0$ as $x\to\infty$.
As you point out, eq. $(1\text a)$ is the differential equation for the parabolic cylinder functions. However, as the DLMF dutifully points out, there is a wide set of possible solutions to this equation. They are of course all equivalent, and any linearly-independent pair serves as a basis, but they are not all equally convenient for all situations.
For this situation, the most convenient solutions are the even and odd solutions given by Wikipedia, which read \begin{align} \psi_E^{\rm(even)}(x) & = e^{-x^2/2}\;_{1}F_{1}\mathopen{}\left(\tfrac14-\tfrac12 E ; \tfrac12 ; x^2 \right)\mathclose{} \\ \psi_E^{\rm(odd)}(x) & = x \, e^{-x^2/2}\;_{1}F_{1}\mathopen{}\left(\tfrac34-\tfrac12 E ; \tfrac32 ; x^2 \right)\mathclose{} \end{align} when adapted to this notation, in terms of the confluent hypergeometric function $_{1}F_{1}$.
Note that these functions almost always blow up at infinity, with the only exceptions being (of course) when $E=\frac12+n$ for an integer $n$ that shares the parity of the eigenfunction (in which case the eigenfunctions obviously reduce to the usual harmonic-oscillator eigenfunctions). However, since we only care about the parabolic-cylinder equation over a finite interval, this is completely irrelevant to our purposes.
This gives us the explicit answer to your question $(1)$, in an explicit (special-function) closed form for the eigenfunctions, which is only complemented by the $|x|>a$ behaviour, $$ \psi_E^{\rm(exp)}(x) = e^{-\sqrt{2(V_0-E)} \, x}. $$
The only thing that's missing is the eigenvalue, which we get by requiring the values of the eigenfunction and its derivative to match at the boundary point, i.e. that \begin{align} A\psi_E^{\rm(even)}(a) & = B \psi_E^{\rm(exp)}(a) \\ A{\psi_E^{\rm(even)}}'(a) & = B {\psi_E^{\rm(exp)}}'(a) \end{align} (and ditto for the odd eigenfunctions). The easiest way to solve this is to realize that it's asking for a nonzero solution $(A,B)$ to the linear system \begin{align} \begin{pmatrix} \psi_E^{\rm(even)}(a) & \psi_E^{\rm(exp)}(a) \\ {\psi_E^{\rm(even)}}'(a) & {\psi_E^{\rm(exp)}}'(a) \end{pmatrix} \begin{pmatrix} A \\ B \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \end{align} and that this can only happen if the system is linearly dependent, which is most cleanly expressed by the vanishing determinant (a.k.a. Wronskian) $$ W_{a,V_0}^{\rm(even)}(E) = \det\mathopen{} \begin{pmatrix} \psi_E^{\rm(even)}(a) & \psi_E^{\rm(exp)}(a) \\ {\psi_E^{\rm(even)}}'(a) & {\psi_E^{\rm(exp)}}'(a) \end{pmatrix}\mathclose{} =0 . \tag 2 $$ For fixed $a$ and $V_0$, this is a function of $E$ that you can easily calculate. Here it is for one particular incarnation, with $a=3$ and $V_0=10$:
Graphically, it's clear that there will be four even bound states for this set of parameters, and the corresponding eigenvalues can be obtained via numerical solution of the transcendental equation $(2)$. There is, of course, no closed-form solution to this transcendental equation in the general case.
A final, minor note: I have no idea what you mean by
That or there are very few or zero bound states, which doesn't match up with my numerical simulation.
You are always guaranteed one bound state, but you are never guaranteed more, and you are never guaranteed an infinite sequence of them. In this particular case, there will indeed be "few" (i.e. a finite number of) bound states. To the extent that they "don't match up" with a numerical simulation, you're essentially guaranteed to be over-interpreting what the numerical simulation tells you. Numerical simulations can never model a true continuum, which means that it takes special care to distinguish bound from continuum states in a numerical eigensolver. However, since you have not explained in depth what you mean by that, it's impossible to say more.