Consider the function $f(y) = (y^2 + c)e^y$. Then $f$ is continuous and has limiting behaviour $\lim_{y \to -\infty} f(y) = 0$ and $\lim_{y \to \infty} f(y) = \infty$. By the intermediate value theorem, there is a solution to $f(y) = x$ for any $x \in (0, \infty)$.
This misses a little bit of information, as if $c < 0$ then there are solutions for some $x < 0$. The same reasoning will work for us. In particular, when $c < 0$ then $f$ takes its minimum at $y_m = -1 + \sqrt{1 - c}$.
We can see this by taking the derivative and setting equal to $0$. The solutions to $f'(y) = 0$ are $y = -1 \pm \sqrt{1 - c}$. The left point $y_M = -1 - \sqrt{1 - c}$ is a maximum. The right point $y_m = -1 + \sqrt{1 - c}$ is a minimum. There are no other extrema. For ease, denote the local minimum $f(y_m) = m$ and the local maximum $f(y_M) = M$. With this, we can classify all solutions.
When $c < 0$, there are solutions to $f(y) = x$ for all $x \geq m$.
- There is exactly one solution to $f(y) = m$.
- For $x \in (m, 0]$ there are exactly two solutions to $f(y) = x$. Further, one solution is in $(y_M, y_m)$ and the other is in $(y_m, \infty)$.
- For $x \in (0, M)$ there are exactly three solutions to $f(y) = x$. One solution is in $(-\infty, y_M)$, another in $(y_M, y_m)$, and the third in $(y_m, \infty)$.
- There are exactly two solutions to $f(y) = M$. One is $y_M$, and the other is in $(y_m, \infty)$.
- There is exactly one solution to $f(y) = x$ for all $x > M$.
One can state analogous analysis for $c \geq 0$.
This has all been for the existence of solutions. Numerically, it is relatively easy to find solutions, especially if we restrict the domain in $y$ to $[y_m, \infty)$. In this region, $f(y)$ is strictly increasing and one can approximate $f(y)$ by $y^2 e^y$, use Lambert's $W$ function, and proceed by bisection to compute solutions.