If you want the states themselves, It's essentialy a question about solving PDEs - you're meant to solve the Schrödinger equation for that potential, i.e.
$$\frac{-\hbar^2}{2m}\nabla^2 \psi + \left[\frac{-A}{y^2+a^2}+\frac{-A}{(y-r)^2+a^2} \right]\psi = E \psi$$
That particular choice of potential doesn't look very solvable analytically though. If you actually wanted the wave functions themselves, you'd need to use some numerical solver.
However, we can still extract useful information - such as "Is there a ground state", and "What is the parity of a given solution". First of all, the variational principle states that for any wavefunction $|\phi\rangle$, $\langle \phi | H | \phi \rangle \ge E_0$, where $E_0$ is the ground state energy. So if you find any old trial state $\phi$ with $\langle \phi | H | \phi \rangle <0$, you have a rigorous bound on the $E_0$ and a proof of the existence of a bound state.
For parity, suppose now that $\psi$ is an exact solution to the above PDE. Then, suppose we invert the parity: Remap $\vec{r} \mapsto -\vec{r}$. Then we end up with
$$\frac{-\hbar^2}{2m}\nabla^2 \psi + \left[\frac{-A}{(-y)^2+a^2}+\frac{-A}{(-y+r)^2+a^2} \right]\psi = E \psi$$
Clearly, there is no change - i.e. any solution of the equation remains a solution under parity change. We can't really say any more about the parity of specific solutions without first finding them though.
