Eigenfunctions and energies of a quadratic-bottom finite square well Consider the potential described by
$$V(x)=\cases{
    \frac{1}{2}m\omega^2x^2 & $\vert x\vert\leq a$\\
    V_0 & otherwise
}$$
where $V_0 > \frac{1}{2}m \omega^2 a^2$, so that we have a potential that looks like this:



*

*Can we write down the solutions to this potential in closed form? 

*Why aren't they the same as in the simple harmonic oscillator case? 


2) is the kicker here. I can write this statement formally in a number of increasingly elaborate ways, but probably the easiest is this: since the generating equation is entirely local (depends only on values of $\psi, \psi^\prime,\psi^{\prime\prime}$, etc. at a single $x$), any solution within the region defined by $\vert x \vert \leq a$ is obeying an equation identical to the simple QHO case. In other words, it sure seems like
$$\left(-\frac{\hbar^2}{2m}\partial^2 + \cases{
    \frac{1}{2}m\omega^2x^2 & $\vert x\vert\leq a$\\
    V_0 & otherwise
}\right) \psi = E\psi$$
should be equivalent to
$$\cases{\left(-\frac{\hbar^2}{2m}\partial^2 + \frac{1}{2}m\omega^2x^2 \right)\psi = E\psi & $\vert x\vert\leq a$\\
\left(-\frac{\hbar^2}{2m}\partial^2 + V_0\right)\psi = E\psi & $\vert x\vert > a$}$$
in other words the equation separates cleanly across the boundary $\vert x\vert=a$. This would in turn would require that $\psi$ be an eigenfunction of each component potential independently within their respective domains, with only the requirement that the eigenvalues match and that the join is continuous and smooth. In particular, that would require $\psi$ match with some QHO solution within the well, $\psi = \cases{\psi^{QHO}_n&$\vert x \vert\leq a$}$ and $E=\left(n+\frac{1}{2}\right)\hbar\omega$ for some eigenfunction of the harmonic oscillator $\psi^{QHO}_n$ with number $n$.
... But this is empirically not true; by which I mean, if you put this into a numerical solver, you will get results that are closer to what we might intuitively expect: the solutions vary smoothly from clean QHO solutions and equal energy intervals for functions bound primarily by the quadratic component, to clean square well solutions at quadratic energy intervals for functions bound primarily by the square well component.

Edit:
If we write the QHO solutions in terms of parabolic cylindrical functions, we can easily write down the second set of normally-diverging solutions: we have $\psi^{QHO}_n(x) \propto D_n(x \sqrt{2})$, and so 
$$\psi = a\, D_n(x\sqrt{2}) + b\, D_{-n-1}(i x\sqrt{2})$$
where for positive integer $n$ the first term generates the solutions we're familiar with. However, in general, the second set of solutions are complex -- and outside the well, we have a decaying exponential $e^{\pm kx}$ that is always pure real and whose derivative is also pure real. So it sure seems like the second set of solutions just gets killed off right away, and we're back where we started...? That or there are very few or zero bound states, which doesn't match up with my numerical simulation.
 A: 
  
*
  
*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.
A: 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.  
