How to solve the 10000th eigenvalue of the anharmonic oscillator? Given a certain Hamiltonian, for example,
$$ H = -\frac{1}{2}\frac{\partial^2}{\partial x^2 } + x^4 .   $$ 
, what methods can we use to approximate the $n$th eigenvalue, for very large $n$?
For small $n$ -- the ground state and the first few low-energy states -- the eigen-energies are easy to calculate, say, with numerical methods, perturbation theory, or even the variational method. But the standard literature of QM doesn't cover the calculation of very high eigen-energies.
Is there any reliable method to approximate $E_n$ for large $n$, preferably with some control over the error?
 A: As @OON mentioned, the WKB is the archetypal tool in the cottage industry studying this system, starting with Bender & Woo 1969 and going on to Cornwall & Tiktopoulos 1993 and on and on... 
To leading order, the energy in WKB is well-known to simply amount to the Bohr-Sommerfeld-Wilson quantization condition. You correctly infer that $\hbar$ is superfluous and absorbable in the variables to 1, but let's keep it here,
$$
\left (-\frac{\hbar^2}{2} \partial_x^2 +x^4-E_n\right ) \psi_n(x)=0,
$$
to connect to 101 yr old formulas, namely, the semiclassical phase integral with turning points at $x=\pm E_n^{1/4}$,
$$
\int_{-E_n^{1/4}}^{E_n^{1/4}}dx \sqrt{2(E_n-x^4)}= (n+1/2)\pi \hbar,
$$ 
for n=0,1,2,3,...,10000, ...  (The 1/2 is the so-called Maslov correction, cf Hall ch 15.2 .)
This is immediately solved to 
$$
E_n=\left (  \frac{(n+1/2)\pi \hbar}{\sqrt{2} C}               \right )^{4/3}  ~,
$$
where $C\equiv \int_{-1} ^1 dy \sqrt{1-y^4}= B(1/4,3/2)/2= \sqrt{\pi}~ \Gamma (5/4)/\Gamma (7/4) \sim 1.748$, so might skip the 1/2 for your chosen n, for which the "many wavelengths in the well" part of the approximation works all too well. 
In fact, if you pursued the references, you'd see that this system is the prototype laboratory for WKB asymptotics methods.
