Apart from a constant normalization, here are the wavefunctions of the QHO.  Equation 2.85 Griffiths.

f[x_,n_] = 1/Sqrt[2^n n!] HermiteH[n,x] Exp[-x^2]

Plot[f[x,0],{x,-3,3},PlotRange->{{-3,3},{-1,1}}]
Plot[f[x,1],{x,-3,3},PlotRange->{{-3,3},{-1,1}}]
Plot[f[x,2],{x,-3,3},PlotRange->{{-3,3},{-1,1}}]
Plot[f[x,3],{x,-3,3},PlotRange->{{-3,3},{-1,1}}]
Plot[f[x,4],{x,-3,3},PlotRange->{{-3,3},{-1,1}}]

And the corresponding probabilities

Plot[f[x,0]^2,{x,-3,3},PlotRange->{{-3,3},{0,1}}]
Plot[f[x,1]^2,{x,-3,3},PlotRange->{{-3,3},{0,1}}]
Plot[f[x,2]^2,{x,-3,3},PlotRange->{{-3,3},{0,1}}]
Plot[f[x,3]^2,{x,-3,3},PlotRange->{{-3,3},{0,1}}]
Plot[f[x,4]^2,{x,-3,3},PlotRange->{{-3,3},{0,1}}]

p0 = Plot[f[x,0]^2,{x,-3,3},PlotRange->{{-3,3},{0,4}}];
p1 =  Plot[f[x,1]^2+1,{x,-3,3},PlotRange->{{-3,3},{0,4}}];
p2 =  Plot[f[x,2]^2+2,{x,-3,3},PlotRange->{{-3,3},{0,4}}];
p3 =  Plot[f[x,3]^2+3,{x,-3,3},PlotRange->{{-3,3},{0,4}}];
pquad = Plot[x^2, {x,-3,3}];
Show[p0,p1,p2,p3,pquad]