(* Mathematica program for Problem 6.10

  date : 23/11/2009
  author: Jorge C. Romão
  email : jorge.romao@ist.utl.pt
*)

(* Question 1 *)
(* Answer given in the problem *)

(* Question 2: Harmonic Oscilator *)

En=Function[n,hbar w (n + 1/2)];

(* <n|x|0> only different from zero for n=1 *)

Test=(En[1]-En[0]) (Sqrt[hbar/(2 m w)])^2;
Print["Question 2"]
Print["Harmonic Oscilator: Sum = ",Test]

(* Question 3: Particle in an Infinite Well *)

En=Function[n,hbar^2 Pi^2/(2 m a^2) n^2];
u=Function[{n,x},Sqrt[2/a] Sin[n Pi x /a]];
xn0=Function[n,Integrate[u[n,x] x u[1,x],{x,0,a}]];
Print["Question 3"]
Print["Infinite Well: <n|x|0> = ",xn0[n] /. {Sin[n Pi]->0,Cos[n Pi]->(-1)^n}]
Print["Infinite Well: Sum = ",Sum[(En[n]-En[1]) xn0[n]^2,{n,1,Infinity}]]

(* Question 4: Half Harmonic Oscillator (Problema 4.12) *)

xinv=Sqrt[m w /hbar];
u=Function[{n,x},Sqrt[xinv] (1/Pi)^(1/4) 1/Sqrt[2^n Factorial[n]] 
	   HermiteH[n,x xinv] E^(-(x xinv)^2/2)];
uhalf=Function[{m,x}, Sqrt[2] u[2 m +1,x]];
xn0 = Function[n, Integrate[uhalf[n, x] x uhalf[0, x], {x, 0, Infinity}, 
    Assumptions -> {m > 0, w > 0, hbar > 0}]];
Print["Question 4"]
Print["Half Harmonic Oscillator: <n|x|0> = ",xn0[n]]
xavgsq=Function[n,(4^n*hbar)/(m*w*(1 + 2*n)!*Gamma[3/2 - n]^2)];
En=Function[n,hbar w (2 n +1 + 1/2)];
g = Function[n, (En[n] - En[0]) xavgsq[n]];
Print["Half Harmonic Oscillator: Sum =",Sum[g[n], {n, 1, Infinity}]]

(* Question 5: Linear Potential (Problem 4.17) *)
Print["Question 5"]

(* The eigenfunctions are

 u_n(y)= C AiryAi[y-eps_n]
 
where 

y=(2 m lbd /hbar^2)^(1/3) x

eps_n = (2 m /(hbar^2 lbd^2))^(1/3) E_n

and eps_n are the zeros of the AiryAi function. These can be obtained 
for instance using *)

ans=FindRoot[AiryAi[x]==0,{x,-3}];
Print["First Energy in units of ",((hbar^2 lbd^2)/( 2m))^(1/3)," = ",-x /. ans]

(* but Mathematica has already a function with the zeros calculated *)

zeros=Table[N[AiryAiZero[i]],{i,500}];
epsn=Table[ -zeros[[n]],{n,1,500}];
Normsq=Table[Quiet[NIntegrate[AiryAi[y-epsn[[i]]]^2,{y,0,Infinity}]],{i,1,500}];
u=Function[{n,y},AiryAi[y-epsn[[n]]]/Sqrt[Normsq[[n]]]];

CheckNormalization:=Table[Quiet[NIntegrate[u[n,y]^2,{y,0,Infinity}]],{n,1,500}]

yn0 = Table[Quiet[NIntegrate[u[n,y] y u[1,y],{y,0,Infinity}]],{n,1,500}];
g= Function[n, (epsn[[n]]-epsn[[1]])  yn0[[n]]^2] ;

Print["Linear Potential: Take 10,20,100,500 terms in the Sum"]
Print["Linear Potential: Sum(10)  = ",hbar^2/(2 m)," ",Sum[g[n], {n, 1, 10}]]
Print["Linear Potential: Sum(20)  = ",hbar^2/(2 m)," ",Sum[g[n], {n, 1, 20}]]
Print["Linear Potential: Sum(100) = ",hbar^2/(2 m)," ",Sum[g[n], {n, 1, 100}]]
Print["Linear Potential: Sum(500) = ",hbar^2/(2 m)," ",Sum[g[n], {n, 1, 500}]]

(* Question 6: Negative Delta Function plus Infinite Well  *)
Print["Question 6"]

(* Initial Values *)
(* Choose lambda in order to exist always a negative energy bound state *)
Clear[lbd]
lbd=3;

(* Width of the well 2 L = 2 a eta  where a is the parameter in the delta 
   function 
 *)

uodd=Function[{n,eta,x},Sqrt[1/eta] Sin[n Pi x /eta]];

(* Energy in units of \hbar^2/(2 m a^2) *)
Eodd=Function[{n,eta},Pi^2 /(4 eta^2) (2 n)^2]  ;

u0=Function[{x,eta}, y0= y /. FindRoot[Tanh[y eta]== 2 y /lbd,{y,1}]; 
	    If[x<-eta,0,If[x<0, Sqrt[(2 y0)/(Sinh[2 y0 eta]- 2 y0 eta)]
			   Sinh[y0 (x +eta)],
	    If[x<eta, Sqrt[(2 y0)/(Sinh[2 y0 eta]- 2 y0 eta)]
	       Sinh[y0 (eta-x)],0]]]];

g=Function[{eta,nmax},y0= y /. FindRoot[Tanh[y eta]== 2 y /lbd,{y,1}];
	   Sum[(Eodd[n,eta]-(-y0^2)) NIntegrate[u0[x,eta] x uodd[n,eta,x], 
           {x,-eta,eta}]^2,{n,1,nmax}]];
	     
g10=Table[{eta,g[eta, 10]},{eta,1,104,5}];
g20=Table[{eta,g[eta, 20]},{eta,1,104,5}];
g50=Table[{eta,g[eta, 50]},{eta,1,104,5}];
g100=Table[{eta,g[eta, 100]},{eta,1,104,5}];

Print["Sum = ",hbar^2/(2m)," g( η)"]
ans=ListPlot[{g10,g20,g50,g100}, AxesLabel -> {"η", "g(η)"}, 
	     PlotStyle->PointSize[0.01], PlotRange->All]