I get a problem with the following nonlinear second order ODE:
$y^2\dfrac{d^2u(y)}{dy^2}+(\Sigma-\rho y^\lambda-\Omega y)y\dfrac{du(y)}{dy}-(\alpha\rho y^\lambda+\beta \Omega y) u(y)=0$,
where $\Sigma, \rho,\Omega,\alpha,\beta,\lambda$ are constants, and $\lambda$ is positive real number.
I want to ask for the exact solution to it. For example, for the special case that $\rho=0$, one of the exact solutions is $M(\beta,\Sigma,\Omega y)$ in the form of Kummer function.
Another pair of representations of the general solution of the $\rho = 0$ specialization, found using Mathematica 11.3.0, which probably used hypergeometric substitution and the contiguous identities to coerce the resulting expressions into "simpler" functions via additional hypergeometric identities. (In other words, it's a black box, but Gosper's algorithm isn't entirely magical.)
$$ u(y) = C_1 \frac{ 1}{\Gamma (\beta ) \Gamma (\beta -\Sigma +1)} G^{2 \, 1}_{1 \, 2}\left( \Omega y \;\middle|\; \begin{matrix} 1-\beta \\ 0, 1 - \Sigma \end{matrix} \right) + C_2 \Gamma (\Sigma -\beta ) G^{2 \, 1}_{2 \, 2}\left( \Omega y \;\middle|\; \begin{matrix} 1-\beta \\ 0, 1 - \Sigma \end{matrix} \right) \text{,} $$ where $C_1$ and $C_2$ are the (arbitrary) constants of integration and $G^{m \, n}_{p \, q}\left(z \;\middle|\; \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \right) $ is the Meijer $G$ function. (Also at MathWorld.)
Using somewhat less exotic functions, $$ u(y) = C_1 U(\beta ,\Sigma , \Omega y ) + C_2 L_{-\beta }^{(\Sigma -1)}( \Omega y ) \text{,} $$ where $C_1$ and $C_2$ are the (arbitrary) constants of integration, $U(a,b,z)$ is the Tricomi (confluent) hypergeometric function (also at MathWorld) and $L^{(\alpha)}_n(x)$ is a generalized Laguerre polynomial (also at MathWorld, without parentheses on the superscript).
There are enough minor variations in the notation for these functions that I might have minor errors in the above. I've combed through these twice, but it's still possible something slipped past me.
Mathematica was unable to solve the version with general $\rho$. I made some progress for rational $\lambda$, where the change of variable $y \mapsto \mathrm{e}^z$ sometimes helps, but that doesn't seem to be relevant given the current version of the Question.
Thoughts on rational $\lambda$, as requested in comments ...
Generally speaking, the $y^\lambda$s produce singular behaviour at $z=0$, which can cause some problems with various solution techniques, especially numerical techniques. So suppose $\lambda = p/q$ for $p,q \in \mathbb{Z}$, $p>0$, $q>0$, and $\gcd(p,q) = 1$. Then a natural choice is to pretend that $y$ and $y^\lambda = y^{p/q}$ are actually just powers of $y^{1/q}$ so that, taking $z = y^{1/q}$, the equation is $$ (q-1)z u''(z) - q(\Sigma - \rho z^p - \Omega z^q) z u'(z) + q^2(\alpha \rho z^p + \beta \Omega z^q)u(z) = 0 \text{.} $$ This isn't suddenly solvable, and still risks singular behavious at $z = 0$. However, it is better behaved when we substitute in a power series for $u$ and, as long as $p$ and $q$ are not too large, the resulting recurrence equations for the coefficients of the series may not be too awful.
A better idea is to push the singularity at $y=0$ infinitely far away. The change of variables $y = \mathrm{e}^z$ does this, pushing the singularity off to the left in the complex plane. Now the equation is $$ u''(z) - (\Sigma - \rho \mathrm{e}^{\lambda z} - \Omega \mathrm{e}^z)u'(z) + (\alpha \rho \mathrm{e}^{\lambda z} + \beta \Omega \mathrm{e}^z) u(z) = 0 \text{.} $$ Power series solutions do not have their radii of convergence limited by whatever is going on at $y = 0$. Again, not suddenly solved, but better behaved for several kinds of attacks.