Exploration of the quantized harmonic oscillator

First look at the normalized solutions of the TISE, 4.12 and 4.13 in the book

c[n_] =  1/Sqrt[Sqrt[Pi] 2^n n!]
ϕ[n_, y_]  =  c[n] Exp[-y^2/2] HermiteH[n, y]

Here the first few Hermite polynomials
Do[
Plot[
            HermiteH[i - 1, x], {x, -1, 1}], 
{i, 5}]

Mathematca can solve for the eigenvalues (energy levels)

n = 0
Solve [
      1/2 {y^2  ϕ[n, y]  -  
                D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, 
              y]  ==  0]

n = 1
Solve [
      1/2 {y^2  ϕ[n, y]  -  
                D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, 
              y]  ==  0]

n = 2
Solve [
      1/2 {y^2  ϕ[n, y]  -  
                D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, 
              y]  ==  0]

n = 3
Solve [
      1/2 {y^2  ϕ[n, y]  -  
                D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, 
              y]  ==  0]

n = 4
Solve [
      1/2 {y^2  ϕ[n, y]  -  
                D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, 
              y]  ==  0]

Do[
Print[ϕ[n, y]], 
{n, 0, 
        5}]

But it's not hard to do it yourself, just plot the trial solution in and simplify.  In each equation below,
if you set the result to zero, the eigenvalue pops right out.  Remember, the result has to be zero for all y.

n = 0;
1/2 {y^2  ϕ[n, y]  -  
              D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, y]  // 
      FullSimplify

n = 1;
1/2 {y^2  ϕ[n, y]  -  
              D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, y]  // 
      FullSimplify

n = 2;
1/2 {y^2  ϕ[n, y]  -  
              D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, y]  // 
      FullSimplify

n = 3;
1/2 {y^2  ϕ[n, y]  -  
              D[ϕ[n, y], {y, 2}]}  -  ϵ ϕ[n, y]  // 
      FullSimplify

Clear[a, n, x]

    ψ[a_, n_, x_]  =  
      1/Sqrt[Sqrt[Pi] 2^n n! a] Exp[-x^2/(2  a^2)] 
        HermiteH[n, x/a]

Here are the first few wavefunctions.

Do[
Plot[ψ[1, i, x], {x, -5, 
          5}], 
{i, 0, 5}]

And here is one of them (n = 
                5) for many snapshots .   To animate select the bracket on the right that encloses all the snapshot then do control - y

i = 0;
nsteps = 25;
n = 5;
e  =  n  +  1/2
Do[ 
      t = i/nsteps; 
Plot[ψ[1, n, x] 
          Cos[e t], {x, -5, 5}, 
        PlotRange -> {{-5, 5}, {-1, 1}}], 
{i, 
        0, nsteps}]