Question: Using the WKB method, provide an approximation for the eigenvalue, $\lambda$, of the problem: $$y'' + \pi^2\lambda y(1+2\cos \pi x)\sin^2(\pi x/2) = 0,~0\le x\le 1,~y(0)=y(1)=0.$$ Compute the first five eigenvalues.
My approach: The ODE is of the form $$Y_{XX}+f(\epsilon X)Y = Y_{XX}+f(x)Y = 0,$$ where $\epsilon = (1/\lambda)^{1/2}$, $X=x/\epsilon$, $Y(X) = Y(x/\epsilon) = y(x),$ and $f(x) = \pi^2(1+2\cos\pi x)\sin^2(\pi x/2).$ Now, $$f(x) \begin{cases}>0,&0<x<2/3,\\ =0,&x=2/3,\\<0,&2/3<x<1.\end{cases}$$ Therefore, we have a turning point at $x=2/3$. As $x\to 2/3$, $$f(x) \sim f'(2/3)(x-2/3)+O\left[(x-2/3)^2\right] \sim \frac{\pi^3}{2}[2\sin(2\pi x)-\sin(\pi x)]_{x=2/3}(x-2/3) = -\frac{3\sqrt{3}\pi^3}{4}(x-2/3)\neq 0.$$ Therefore, $x=2/3$ is a first order turning point.
Now, WLOG, assuming $\lambda > 0$ and $g(x) \equiv \lambda f(x)$, we have $$g(x) \begin{cases}>0,&0<x<2/3,\\ =0,&x=2/3,\\<0,&2/3<x<1.\end{cases}$$ Therefore, by WKB, when $0<x<2/3$, we have the oscillatory solution: $$y(x;\epsilon) \sim \frac{1}{g(x)^{1/4}}(a\cos \theta + b\sin\theta)=\frac{1}{\left[\lambda\pi^2\sin^2(\pi x/2)(1+2\cos\pi x)\right]^{1/4}}(a\cos\theta + b\sin \theta),$$ where $$\theta(x) = \pi|\lambda|^{1/2}\left|\int_{\frac{2}{3}}^x[1+2\cos (\pi x')]^{1/2}\sin(\pi x'/2)dx'\right|.$$ Now, using numerical integration $\theta(0) \approx 0.3417\pi|\lambda|^{1/2}.$
Similarly, by WKB, when $2/3<x<1$, we have the exponential solution: $$y(x;\epsilon) \sim \frac{1}{\sqrt{2}[-g(x)]^{1/4}}\left[(a-b) e^{\Phi} + 2(a+b)e^{-\Phi}\right]=\frac{1}{\sqrt{2}|\lambda|^{1/4}\left[\pi^2\sin^2(\pi x/2)(1+2\cos\pi x)\right]^{1/4}}\left[(a-b) e^{\Phi} + 2(a+b)e^{-\Phi}\right],$$ where $$\Phi(x) = \pi|\lambda|^{1/2}\left|\int_{\frac{2}{3}}^x|1+2\cos (\pi x')|^{1/2}\sin(\pi x'/2)dx'\right|.$$ Now, using numerical integration $\Phi(1) = 0.25\pi |\lambda|^{1/2}$. This implies $$\frac{1}{\sqrt{2}\pi^{1/2}|\lambda|^{1/4}}[(a-b)e^{\Phi(1)} + 2(a+b)e^{-\Phi(1)}]=0\implies (a-b)e^{2\Phi(1)} + 2(a+b)=0$$ $$\implies b = a\left(\frac{e^{2\Phi(1)}+2}{e^{2\Phi(1)}-2}\right).$$
At this point, notice that it is not possible to satisfy the boundary condition $y(0) = 0$ for any combination of $a$ and $b$ in the oscillatory solution, since $g(0)=0$. Therefore, we have a boundary layer near $x=0$.
Now, using Taylor expansion, near $x = 0$, $g(x)\sim \lambda \pi^2 x^2 + O(x^4)$. Therefore, in the boundary layer, the ODE can be rewritten as $$w'' + \lambda \pi^2 x^2 w = w'' + (|\lambda|^{1/2} \pi)^2 x^{4-2} w = 0,$$ whose solution is $$w(x) = \sqrt{x}\mathcal{C}_{1/4}\left(\frac{1}{2}|\lambda|^{1/2} \pi x^2\right)\sim \frac{2}{\pi|\lambda|^{1/4}x^{1/2}}\cos\left(\frac{1}{2}|\lambda|^{1/2} \pi x^2 - \frac{3\pi}{8}\right),$$ where $\mathcal{C}_{\nu}$ is a Bessel function of order $\nu$.
I am not exactly sure how to proceed after this. Can someone please help me? Thanks in advance.
Edit: Alternatively, note that $y(0)=0$ only when $$a \cos[\theta(0)] + b\sin[\theta(0)] = 0$$ $$\implies a\cos[\theta(0)] + a\left(\frac{e^{2\Phi(1)}+2}{e^{2\Phi(1)}-2}\right)\sin[\theta(0)] = 0$$ $$\implies \cos[\theta(0)] + \left(\frac{e^{2\Phi(1)}+2}{e^{2\Phi(1)}-2}\right)\sin[\theta(0)]=0,$$ since $a=0$ implies $b=0$, which would lead to the trivial solution. Carrying on $$\tan\left(0.3417 \pi|\lambda|^{1/2}\right) = \frac{2-e^{0.5\pi|\lambda|^{1/2}}}{2+e^{0.5\pi|\lambda|^{1/2}}}\sim -1$$ $$\implies|\lambda_n| \sim 8.569\left(n-\frac{1}{4}\right)^2,~n\in\mathbb{Z_{\ge 0}}.$$