Harmonic Oscillator

Show[Plot[(x^2),{x,-3,3},PlotRange->{{-3,3},{0,5}}],ListLinePlot[{{{-(Sqrt[0.5]),0.5},{(Sqrt[0.5]),0.5}},
{{-(Sqrt[1.5]),1.5},{(Sqrt[1.5]),1.5}},
{{-(Sqrt[2.5]),2.5},{(Sqrt[2.5]),2.5}},
{{-(Sqrt[3.5]),3.5},{(Sqrt[3.5]),3.5}},
{{-(Sqrt[4.5]),4.5},{(Sqrt[4.5]),4.5}}}]]

eqn=-ℏ^2/(2μ)ψ''[x]+k x^2/2 ψ[x]==Ev ψ[x]

Ev:=(-ℏ^2/(2μ)ψ''[x]+k x^2/2 ψ[x])/ψ[x]

α=Sqrt[k μ]/ℏ

ψ[x_]=Exp[-α x^2 / 2]

Ev

Simplify[%,{k > 0, μ > 0, ℏ > 0}]

PowerExpand[%]

ψ[x_]=x Exp[-α x^2 / 2]

Simplify[Ev,{k > 0, μ > 0, μ > 0}]

PowerExpand[%]

Cnorm[v_]:=Sqrt[(1/(2^v v!))](α/Pi)^(1/4) 
Cnormv[v_]:=Sqrt[(1/(2^v v!))]

Cnorm[0]

Cnormv[0]

Herm[0,y_]=1

Herm[1,y_]:= 2y

Herm[v_,y_]:=2 y Herm[v-1,y]-2(v-1)Herm[v-2,y]

Herm[3,y]

Herm[2,y]

ψv[v_,y_]:=Cnorm[v] Herm[v,y] Exp[-y^2/2]

ψv[0,y]

ψv[1,y]

Integrate[ψv[0,y]^2 / Sqrt[α],{y,-Infinity,Infinity}]

Integrate[ψv[1,y]^2 / Sqrt[α],{y,-Infinity,Infinity}]

Integrate[ψv[2,y]^2 / Sqrt[α],{y,-Infinity,Infinity}]

Integrate[ψv[3,y]^2 / Sqrt[α],{y,-Infinity,Infinity}]

Plot[{Cnormv[0]Herm[0,y] Exp[-y^2/2],Cnormv[1]Herm[1,y] Exp[-y^2/2],Cnormv[3]Herm[3,y] Exp[-y^2/2],Cnormv[10]Herm[10,y] Exp[-y^2/2]},{y,-3,3}]

Manipulate[
Plot[{Cnormv[n]Herm[n,y]Exp[-y^2/2],(((Cnormv[n]Herm[n,y]Exp[-y^2/2]))^2)},{y,-6,6},AxesLabel->{"y","ψ(y) and ψ(y)^2"},PlotLabel -> "Harmonic Wavefunctions",  PlotRange->{{-6,6},{-1,1}},PlotLegends->{"ψ","(ψ^2)"}],
{n,0,12,1}]

Evib[v_]:=(-ℏ^2/(2μ)D[ψv[v,Sqrt[α]x],{x,2}]+k x^2/2 ψv[v,Sqrt[α]x])/ψv[v,Sqrt[α]x]/(ℏ Sqrt[k/μ])

condition= {k > 0, μ > 0, ℏ > 0}

Do[Print["v = ",n,"  E = ",Simplify[Evib[n],condition]," hν"],{n,0,3}]

avex=Integrate[y ψv[0,y]^2 / Sqrt[α],{y,-Infinity,Infinity}] / Sqrt[α]

avex2=Integrate[y^2 ψv[0,y]^2 / Sqrt[α],{y,-Infinity,Infinity}] / α

△x=Sqrt[avex2-avex^2]

avep=-ℏ Integrate[ ψv[0,y] D[ψv[0,y],y] ,{y,-Infinity,Infinity}]

avep2=-ℏ^2 Integrate[ ψv[0,y] D[ψv[0,y],{y,2}]  Sqrt[α],{y,-Infinity,Infinity}]

△p=Sqrt[avep2-avep^2]

△x △p

Simplify[△x △p,condition]

eqny=1/2 k (y/Sqrt[α])^2==(v + 1/2)ℏ Sqrt[k/μ]

Solve[eqny,y]

Simplify[%,condition]

ytp[v_]:= Sqrt[1+ 2 v]

Pforb[v_]:=2 Integrate[ψv[v,y]^2 ,{y,ytp[v],Infinity}] /Sqrt[α]

Pforb[0]

N[%]

Pforb[1]

N[%]

N[Pforb[2]]

Pforb[1]

N[Pforb[3]]

Pforb[3]

Pforb[4]

N[%]

Do[Pf[n]=N[Pforb[n]],{n,0,6}]

DiscretePlot[Pf[n],{n,0,6},PlotRange->{{0,6},{0,0.2}}]

Show[Plot[(((1-(E^(-x))))^2),{x,-0.8,5},PlotRange->{{-1,5},{0,1.2}}],ListLinePlot[{{{-Log[1+(Sqrt[0.19])],0.19},{-Log[1-(Sqrt[0.19])],0.19}},
{{-Log[1+(Sqrt[0.51])],0.51},{-Log[1-(Sqrt[0.51])],0.51}},
{{-Log[1+(Sqrt[0.75])],0.75},{-Log[1-(Sqrt[0.75])],0.75}},
{{-Log[1+(Sqrt[0.91])],0.91},{-Log[1-(Sqrt[0.91])],0.91}},
{{-Log[1+(Sqrt[0.99])],0.99},{-Log[1-(Sqrt[0.99])],0.99}}}]]

Evibm[v_]:=(v+(1/2))hν-(((v+(1/2)))^2)hν (Subscript[χ, e])

Do[Print["v = ",n,"  E = ",Evibm[n]],{n,0,3}]

hν=12.59; (Subscript[χ, e])=0.0275;

Do[Print["v = ",n,"  E = ",Evibm[n], " kcal/mol"],{n,0,15}]

De=hν / (4 (Subscript[χ, e]))

k=520

a=(Sqrt[k/(2*De*1000*4.184/6.022*^23)])*1*^-10

De(((1-(E^(-a 1))))^2)

Plot[De(((1-(E^(-a x))))^2),{x,-0.8,5},PlotRange->{{-1,3},{0,De}}]

xtm[v_,1]:=((-1)/a)Log[1+(Sqrt[((Evibm[v])/De)])]

xtm[v_,2]:=((-1)/a)Log[1-(Sqrt[((Evibm[v])/De)])]

ListLinePlot[Table[{{xtm[n,1],Evibm[n]},{xtm[n,2],Evibm[n]}},{n,0,5}]]

Show[Plot[De(((1-(E^(-a x))))^2),{x,-0.8,5},PlotRange->{{-1,3},{0,De+10}},PlotLabel -> "Vibrational Energy Level of Hydrogen Molecule",AxesLabel->{Angstrom,kcal/mol}],ListLinePlot[Table[{{xtm[n,1],Evibm[n]},{xtm[n,2],Evibm[n]}},{n,0,15}]],
ListLinePlot[{{xtm[20,1],De},{3,De}},PlotStyle->Dashed]]