% WKB approximations to eigensolutions of the harmonic oscillator
%   y" + (lambda-x^2)y= 0, y->0 as |x|->infty
%   lambda_n = 2*n+1

  n = 6;   % Index of eigenvalue/eigenfn
  x = -5:0.1:5;
  
  % Calculate Hermite polynomial Hn(x) of order n using recursion
  %  Hk(x) = 2*x*Hkm1(x) - 2*km1*Hkm2(x)
  
  H0 = 1;
  H1 = 2*x;
  switch n
   case 0
    Hn = H0;
   case 1
    Hn = H1;
   otherwise
    Hkm2 = H0;
    Hkm1 = H1;
    km1 = 1;
    for k = 2:n
      Hk = 2.*x.*Hkm1 - 2*km1*Hkm2;
      Hkm2 = Hkm1;
      Hkm1 = Hk;
      km1 = k;
    end
    Hn = Hk;
  end
  
  yex = Hn.*exp(-0.5*x.^2);
  yex = yex/yex(x==0);  % Normalize to 1 at x=0 (requires n even)
  
  % WKB solution
  %  y = A*m^(-1/4) cos(phi - pi/4)
 
  lambda = 2*n+1;
  delta = 0.1;
  xosc = x(abs(x)<sqrt(lambda)-delta);
  xi = acos(xosc/sqrt(lambda));
  phi = lambda*(0.5*(pi-xi)+0.25*sin(2*xi));
  m = sqrt(lambda - xosc.^2);
  yWKB = cos(phi - pi/4)./sqrt(m);
  yWKB = yWKB/yWKB(xosc==0);

  subplot(2,2,1)
  plot(x,yex,'b-',xosc,yWKB,'m--')
  xlabel('x')
  ylabel('y')
  ymax = 2;
  ylim([-ymax ymax])
  legend('Exact','WKB')
  title(['Eigenfn of y"+(2n+1-x^2)y=0, n=' num2str(n)])
  hold on
  plot([-sqrt(lambda) -sqrt(lambda)],[-ymax ymax],'r-.') % Plot turning points
  plot([ sqrt(lambda)  sqrt(lambda)],[-ymax ymax],'r-.')
  hold off