Beyond WKB approximation for energies
In the following, I'll describe to you how in principle one can compute higher order corrections to the Bohr-Sommerfeld condition. In order to find higher order corrections to the Bohr-Sommerfeld formula, we need to include higher order corrections of the wavefunction of the form:
$\Psi(x) = \sum_n \hbar^n a_n(x) exp(\frac{i}{\hbar} \int^x \sqrt{2m(E-V(y)) }dy)$
Inserting in the Schrodinger equation and requiring fullfilment to each order of $\hbar$, we obtain:
$a_0 = \frac{Const.}{(2m(E-V(y)))^{\frac{1}{4}} }$
and
$ a_k \Delta S + 2 \nabla a_k. \nabla S = i \Delta a_{k-1}$
where:$ S= \int^x \sqrt{2m(E-V(y)) } dy$. In one dimensions these equations can solved in closed form.
The quantization condition is achieved through the requirement that the wave function phase does not change over a closed loop:
$\oint^x \frac{1}{\Psi} \frac{d\Psi}{dx} dx = 2\pi i n\hbar $
Please see, the following article by A. Voros where this procedure was applied to the anharmonic oscillator. Equation (4.4), describes an explicit solution of the wave function to the first higher order.
Remark: Voros uses the Bargmann representation of the phase space instead of the usual coordinate-momentum representation.