How to qualitatively draw the stationary states of this problem (of course, without solving it in detail)
Start with the obvious facts:
- The potential is symmetric: $V(-x)=V(x)$.
- The potential is bounded and is a potential barrier of finite width.
- The solution for $|x|>x_0$ is a plane wave.
The first fact lets you classify the solutions in two kinds: even and odd ones. This reduces the problem on $x\in\mathbb R$ to two problems on $x\in[0,\infty)$. The even solutions will have homogeneous Neumann boundary condition at $x=0$, while the odd ones will have the Dirichlet boundary there.
The second fact lets you immediately say that there'll be continuous spectrum. Moreover, the third fact tells you that the problem is a free particle scattering problem, where an incident particle can have any energy, so there are no gaps in the spectrum.
From this follow several ways to proceed.
Purely qualitative approach
If you really don't want to solve the problem quantitatively, and only want to qualitatively sketch the eigenstates, without much regard to exact energies and other numeric details, then you can approximate the finite parabolic potential with a square barrier of some width:

If you are able to sketch eigenstates for this potential, you're done.
Quantitative approach
This is actually quite easy if you can program a solver for ODE Cauchy problems. A method even as simple as Euler's should be enough to get started here.
In particular, you use the insight given by the above mentioned obvious facts and start with a Cauchy problem for each parity of the solution you're looking for. For the odd solution you can impose e.g. the condition
$$\psi_{\text{odd}}(0)=0,$$
$$\left.\frac{d\psi_{\text{odd}}(x)}{dx}\right|_{x=0}=1,$$
and for the even solution you can impose
$$\psi_{\text{even}}(0)=1,$$
$$\left.\frac{d\psi_{\text{even}}(x)}{dx}\right|_{x=0}=0.$$
Solving the Cauchy problems with these initial conditions won't give you normalized (in any sense) solutions, but I guess you don't require them. If you do, then note that these states can't be normalized to probability of 1, since they correspond to unbounded motion, so a different normalization scheme is required.
After you do solve the Cauchy problems and extend the odd and even solutions calculated e.g. for $x>0$ to the space of $x<0$, you'll get something like the following solution for the case of $m=1$, $\omega=1$, $x_0=2$.

In particular, how different will the stationary state in this case be from those of the actual harmonic oscillator potential $V(x)=\frac{1}{2}m\omega^2x^2, ~\forall x$?
For this, I guess, the best answer would be to simply do the numerical calculation and plot the solutions side-by-side. It's very fortunate that one of the problems has continuous spectrum, since we can now take any eigenstate of a SHO and find the solution to the modified problem with the same energy. Below is a set of comparisons of the first several eigenstates of SHO (blue) and of the modified equation (orange) with the same values of $m$, $\omega$ and $x_0$ as taken above.
