If you want a quick and dirty way to find the asymptotics of $y(x)$ as $x \to \infty$, you can use the WKB ansatz $y(x) = e^{S(x)}$, with the hypothesis that $S^{(k)}/S' \to 0$ as $x\to \infty$ for all $k>1$. Substituting this form into your equation and keeping only the leading terms at infinity, you find
$$
x^{2n} (S')^{2n} + x^2 = 0 .
$$
The solutions are $S'(x) = (-x^{2-2n})^{1/2n} = (-1)^{1/2n} x^{1/n-1}$ and hence the expected asymptotics for the $2n$ independent solutions are $y_k(x) \sim \exp(e^{i\frac{\pi}{2n}(1+2k)} n x^{1/n})$, $k=0,\ldots,2n-1$. In the case $n=1$, you get the asymptotics $y(x) \sim e^{\pm i x}$, and not $e^{\pm x}$ as your question suggested. You get the $e^{\pm x}$ asymptotics by flipping the sign of the $x^2 y$ term in the equation.
If you want to be more systematic, it helps to apply some transformations to the equation first. If you introduce $z = x^n y^{(n)}$, then the equation is equivalent to
$$
\begin{bmatrix}
x^n y^{(n)} \\
(x^n z)^{(n)}
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
-x^2+\lambda & 0
\end{bmatrix}
\begin{bmatrix}
y \\
z
\end{bmatrix} .
$$
Next, it is convenient to rescale the dependent variables such that the leading power of $x$ on the right-hand side has a diagonalizable matrix as coefficient (if it's actually possible). Here, this looks like
$$
\begin{bmatrix}
x^n (x y)^{(n)} + P(x\frac{d}{dx}) y \\
(x^n z)^{(n)}
\end{bmatrix}
= x
\begin{bmatrix}
0 & 1 \\
-1+\frac{\lambda}{x} & 0
\end{bmatrix}
\begin{bmatrix}
x y \\
z
\end{bmatrix} ,
$$
where now $P(x\frac{d}{dx})$ is a constant coefficient polynomial in its argument, of order less than $n$. The last step is to substitute $x=t^{n}$, where the power of $t$ is chosen so that the resulting equation has asymptotically constant coefficients as $t\to \infty$. Since $x\frac{d}{dx} = \frac{1}{n} t\frac{d}{dt}$, the result is
$$
\begin{bmatrix}
(t^n y)^{(n)} + t^{-n} Q(t\frac{d}{dt}) y \\
z^{(n)} + t^{-n} R(t\frac{d}{dt}) z
\end{bmatrix}
= n^n
\begin{bmatrix}
0 & 1 \\
-1+\frac{\lambda}{t^n} & 0
\end{bmatrix}
\begin{bmatrix}
t^n y \\
z
\end{bmatrix} ,
$$
where again the constant coefficient polynomials $Q$ and $R$ are of order less than $n$. From here the theory of asymptotics of ODEs at irregular singular points tells us that the independent solutions of your equation will have the asymptotic expansions $y_k(t) = t^{-n} Y_k(t) (1 + O(1/t))$, where the last factor stands for an asyptotic series in inverse powers of $t$ and $Y_k(t)$ solves the asymptotic constant coefficient equation. In this case corresponds to $Y_k(t) = \exp(\alpha_k t)$, where the $\alpha_k$ are the different $n$-th roots of the eigenvalues of matrix $n^n [\begin{smallmatrix} 0 & 1 \\ -1 & 0 \end{smallmatrix}]$, namely $\alpha_k = e^{i\frac{\pi}{2n}(1+2k)} n$, $k=0,\ldots,2n-1$, which reproduces the form of the $S(x)$ phase function from the WKB solution. But now we have more detailed information about the asymptotics.
A good reference for the theory of asymptotics of ODEs with meromorphic coefficients is the following book. However, it deals with higher order equations by reducing them to first order systems, so the details of the calculational steps that I explained above will be slightly different there.
Wasow, W., Asymptotic expansions for ordinary differential equations, New York: Dover. (1965). ZBL0169.10903.
