WKB approximation to the harmonic oscillator wave functions

We take m=1, ω=1, ℏ=1 and work out WKB wave function

First, classically allowed region

Integrate[(Sqrt[(2  n + 1) - (Q^2)]), Q]

Limit[%, Q -> -(Sqrt[2  n + 1])]

u = (1/(((2  n + 1 - (Q^2)))^(1/4))) Cos[%% - % - (π/4)]

Next, classically forbidden region

Integrate[(Sqrt[-(2  n + 1) + (Q^2)]), Q]

Limit[%, Q -> -(Sqrt[2  n + 1])]

v1 = 
      (1/2) (1/(((-(2  n + 1) + (Q^2)))^(1/4))) 
        Exp[-% + %%]

Integrate[(Sqrt[-(2  n + 1) + (Q^2)]), Q]

Limit[%, Q -> (Sqrt[2  n + 1])]

v2 = 
      (1/2) (1/(((-(2  n + 1) + (Q^2)))^(1/4))) 
        Exp[% - %%]

n=20

uplot = 
      Plot[u /. {n -> 20}, {Q, -(Sqrt[41]), (Sqrt[41])}, PlotPoints -> 200]

v1plot = 
      Plot[v1 /. {n -> 20}, {Q, -2 (Sqrt[41]), (Sqrt[41])}, 
        PlotPoints -> 200]

v2plot = 
      Plot[v2 /. {n -> 20}, {Q, (Sqrt[41]), 2 (Sqrt[ 41])}, PlotPoints -> 200]

a1plot = 
      Plot[((((π/(((2 (Sqrt[ 2  n + 1])))^(1/3)))))^(1/2)) 
            AiryAi[(((2 (Sqrt[ 2  n + 1])))^(1/3)) 
                (-Q - (Sqrt[2  n + 1]))] /. {n -> 20}, {Q, 
          ((-3)/2) (Sqrt[41]), ((-3)/4) (Sqrt[41])}]

a2plot = 
      Plot[((((π/(((2 (Sqrt[ 2  n + 1])))^(1/3)))))^(1/2)) 
            AiryAi[(((2 (Sqrt[ 2  n + 1])))^(1/3)) 
                (Q - (Sqrt[2  n + 1]))] /. {n -> 20}, {Q, 
          (3/4) (Sqrt[41]), (3/2) (Sqrt[41])}]

WKBplot = Show[uplot, v1plot, a1plot, v2plot, a2plot]

N[(u
            /(((((Sqrt[π]) (2^n) n!))^(-1/2)) 
              HermiteH[n, Q] (E^(-(Q^2)/2)))) /. {n -> 20} /. {
          Q -> 0}]

trueplot = 
      Plot[1.253407199896294` 
            ((((Sqrt[π]) (2^n) n!))^(-1/2)) 
            HermiteH[n, x] (E^(-(x^2)/2)) /. {n -> 20}, {x, 
          -2 (Sqrt[41]), 2 (Sqrt[ 41])}, PlotPoints -> 200]

Show[WKBplot, trueplot]

n=5

uplot = 
      Plot[u /. {n -> 5}, {Q, -(Sqrt[11]), (Sqrt[11])}, PlotPoints -> 200]

v1plot = 
      Plot[v1 /. {n -> 5}, {Q, -2 (Sqrt[11]), (Sqrt[11])}, 
        PlotPoints -> 200]

v2plot = 
      Plot[-v2 /. {n -> 5}, {Q, (Sqrt[11]), 2 (Sqrt[ 11])}, 
        PlotPoints -> 200]

a1plot = 
      Plot[((((π/(((2 (Sqrt[ 2  n + 1])))^(1/3)))))^(1/2)) 
            AiryAi[(((2 (Sqrt[ 2  n + 1])))^(1/3)) 
                (-Q - (Sqrt[2  n + 1]))] /. {n -> 5}, {Q, 
          ((-3)/2) (Sqrt[11]), ((-3)/4) (Sqrt[11])}]

a2plot = 
      Plot[-((((π/(((2 (Sqrt[ 2  n + 1])))^(1/3)))))^(1/2)) 
            AiryAi[(((2 (Sqrt[ 2  n + 1])))^(1/3)) 
                (Q - (Sqrt[2  n + 1]))] /. {n -> 5}, {Q, 
          (3/4) (Sqrt[11]), (3/2) (Sqrt[11])}]

WKBplot = Show[uplot, v1plot, a1plot, v2plot, a2plot]

N[(u
            /(((((Sqrt[π]) (2^n) n!))^(-1/2)) 
              HermiteH[n, Q] (E^(-(Q^2)/2)))) /. {n -> 5} /. {
          Q -> 0.5}]

trueplot = 
      Plot[-1.2550763185968383` 
            ((((Sqrt[π]) (2^n) n!))^(-1/2)) 
            HermiteH[n, x] (E^(-(x^2)/2)) /. {n -> 5}, {x, 
          -2 (Sqrt[11]), 2 (Sqrt[ 11])}, PlotPoints -> 200]

Show[WKBplot, trueplot]