Consider a time indipendent Schrodinger problem: $$\hat{H}\psi_E(p) = E \psi_E(p)$$ with suitable boundary conditions. We know that $\psi_E$ are the eigenfunctions of $\hat{H}$. If we now consider the operator $\hat{H}^2$ we know that the eigenfunctions of $\hat{H}$ are also eigenfunctions of $\hat{H}^2$ but in general the reverse is not true.
Consider the problem $$\hat{H}^2 \psi_\epsilon (p) = \epsilon \psi_\epsilon (p) \, .$$ What is the relation between the two solutions? My intuitive idea is that the eigenfunctions of $\hat{H}$ are eigenfunctions of $\hat{H}^2$ whenever $E=+\sqrt{\epsilon}$, so that starting from the "squared problem" we could, in principle, find the eigenfunctions and eigenvalues of the "non-squared problem" simply restricting the solutions to $$E = +\sqrt{\epsilon}$$ and $$\psi_E(p) = \psi_{\sqrt{\epsilon}}(p) \, .$$
Is this legit or are there inconsistencies when starting from the "squared problem"? Should I also consider the solutions with $E=-\sqrt{\epsilon}$?