Imagine a particle obeying Schrodinger's equation with the following modified harmonic oscillator potential: $$ V(x)=\frac{1}{2}m\omega^2x^2\pm mgx \ \ \ \ \ \ \text{for} \ \ x\ge0 \\ V(x)=\infty \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{for} \ \ x<0 $$
The time-independent Schrodinger equation is then $$ -\frac{\hbar^2}{2m}\frac{d^2\Psi}{dx^2}+\frac{1}{2}m\omega^2x^2\Psi\pm mgx\Psi=E\Psi \ \ \ \ \ \ \text{for} \ \ x\ge0 $$ which is equivalent to $$ \frac{d^2\Psi}{dz^2}-(z^2-K)\Psi=0 \ \ \ \ \ \ \ \ \ \text{for} \ \ z\ge \pm L $$ where $$ z=\sqrt{\frac{m\omega}{\hbar}}(x\pm\frac{g}{\omega^2}) \\ K=\frac{2E}{\hbar \omega}+\frac{mg^2}{\hbar\omega^3} \\ L=\sqrt{\frac{mg^2}{\hbar\omega^3}} $$
The essential boundary conditions for this problem are $$ \Psi(z\le\pm L)=\Psi(\infty)=0 $$ Based on the analysis I have done so far, it seems as though the right-hand boundary condition (at infinity) should impose the usual normalization constraints on the solution space. Therefore, $$ K_n=2n+1 \rightarrow E_n=\frac{\hbar\omega}{2}(2n+1-\frac{mg^2}{\hbar\omega^3}) \\ \\ \\ \Psi_n(z)=A_nH_n(z)e^{-\frac{z^2}{2}} $$ where $\Psi_n(z)$ is the $n$th eigenfunction of the problem, $A_n$ is some normalization constant for the $n$th eigenfunction, and $H_n(z)$ is the $n$th Hermite polynomial, defined as usual.
However, this solution appears to not have enough degrees of freedom to accommodate the left-hand boundary condition.
Unlike the solutions to the usual half-harmonic oscillator, we can't merely throw out the even-parity solutions. The Hermite polynomials have roots which are determined by their definitions and which cannot be changed (as far as I'm aware) by the parameters of the problem.
This leaves me in (I believe) one of three possible positions:
(a) The eigenfunctions of this problem exsist but are different in kind to the usual harmonic and half-harmonic oscillator solutions. (i.e. I made an analysis mistake.)
(b) The computed eigenfunctions are correct, and there is a way to guarantee the left-hand boundary condition. (i.e. I'm missing a degree of freedom.)
(c) There are no eigenfunctions for this system. Therefore, the system has no states of definite energy, but the left-hand boundary condition can be guaranteed through some convergent series of otherwise valid solutions. (i.e. Something complicated I don't understand is going on here.)
Question: What are the solutions to the above differential equation under the stated boundary conditions?
Short Answer: The eigenfunctions take the form of $D_\nu(x)$, the parabolic cylinder function.
Long Answer: You made an analysis mistake; (a) is correct.
This is a well-stated Sturm-Liouville Problem. Usually, to guarantee the existence of a discrete eigenvalue spectrum, the boundary conditions must apply to a closed interval $[a,b]$ (e.g. where neither $a$ nor $b$ is $\infty$). However, it is a known result that one-dimensional Schrodinger operators with potentials diverging towards infinity on both the left and the right have a discrete eigenvalue spectrum. This means that definite energy states exist which satisfy both boundary conditions.
Instead of making the substitution given in the question, do the following: $$ z=\sqrt{\frac{2m\omega}{\hbar}}x + L \\ \nu=\frac{E}{\hbar\omega}+\frac{L^2}{4}-\frac{1}{2} \\ L=\pm\sqrt{\frac{2mg^2}{\hbar\omega^3}} $$ This transforms the problem into Weber's Equation with the following boundary conditions $$ \frac{d^2\Psi}{dz^2}+(\nu+\frac{1}{2}-\frac{1}{4}z^2)\Psi=0 \ \ \ ; \ \Psi(z=L)=\Psi(z=\infty)=0 $$
This equation has a solution of the form $D_\nu(z)$, the parabolic cylinder function, where only specific (i.e. discrete) values of $\nu$ satisfy the stated boundary conditions.
For example, it is known that $D_n(x)$ can be written in terms of the Hermite polynomials for $n\in \{0,1,2,3,...\}$. Therefore, $L=0$ implies that $E_\nu=(\nu+\frac{1}{2})\hbar\omega$ and that this solution reproduces the familiar odd-parity eigenfunctions cited in the question as solutions to the traditional quantum half-harmonic oscillator.
However, for other values of $L$ simple and exact solutions are not possible. Numerical approximations to within some prescribed tolerance must must be obtained. I used Mathematica to derive the following results.
For the sake of illustration, let's specify $L=-\sqrt{2}$ so that the energy takes on the simple form $E_\nu=\nu\hbar\omega$. In this case, the groundstate of the system corresponds to $\nu\approx 0.234234$, its square-integral gives $\approx 1.66278$, and thus its normalized eigenfunction looks like this: $$ \Psi_{groundstate}(z)=\frac{1}{\sqrt{1.66278}}D_{0.234234}(z) $$
Here are all of the eigenvalues less than $10$ and their corresponding normalization constants for $L=-\sqrt{2}$: $$ \nu\approx 0.234234 \ \ , \ \ N^2\approx 1.66278 \\ \nu\approx 1.69746 \ \ , \ \ N^2\approx 2.51615 \\ \nu\approx 3.28019 \ \ , \ \ N^2\approx 13.3306 \\ \nu\approx 4.92981 \ \ , \ \ N^2\approx 159.493 \\ \nu\approx 6.62244 \ \ , \ \ N^2\approx 3486.3 \\ \nu\approx 8.34549 \ \ , \ \ N^2\approx 122853. \\ $$
Some additional notes:
(i) You may be surprised that the groundstate eigenvalue for this system leads to a ground state energy of less than $\frac{1}{2}\hbar\omega$, the groundstate energy of the unrestricted harmonic oscillator. Restricting the domain ought to increase the groundstate energy, and in a sense, it does. Keep in mind that the way we changed this problem was by adding a linear potential which had the effect of lowering the vertex of the parabola and shifting it to the right. This means that while the absolute energy of the ground state may be lower, its difference with respect to the potential minimum is greater. The key is that $0<0.234234<1$, where $0$ is the eigenvalue for the harmonic oscillator ground state and $1$ is the eigenvalue for the half-harmonic oscillator ground state, so everything is copacetic.
(ii) In a similar vein, if $L>0$, then $\nu>1$ and the groundstate has a higher energy (as expected).
(iii) A related problem exists in which the infinite potential barrier is on the right instead of the left. By symmetry, the solution to this problem is simply $D_\nu(-z)$.
(iv) $D_\nu(z)$ diverges to $\pm\infty$ on the left when $\nu\notin\{0,1,2,3,...\}$, which explains why the solution with arbitrary $\nu$ is not suitable in the unrestricted case.