#----------------------------#
#    Numerischer Test der    #
#    Semiklassischen En's    #
#----------------------------#

#  R-Code  #


beta = 0.01
N = 1000

h0 = diag( (0:N)+1/2 )
h1 = matrix(0,N+1,N+1)

for(n in 0:N)
{
  h1[1+n,1+n] = 6*n*(n+1)+3
  
  if(n+2<=N)
  {
  h1[1+n,1+n+2] = (4*n+6)*sqrt((n+1)*(n+2)) 
  h1[1+n+2,1+n] = h1[1+n,1+n+2]
  }
  if(n+4<=N)
  {
  h1[1+n,1+n+4] = sqrt((n+1)*(n+2)*(n+3)*(n+4))
  h1[1+n+4,1+n] = h1[1+n,1+n+4]
  }
}

h = h0 + beta/16 * h1

res = rev(eigen(h, symmetric=TRUE, only.values=TRUE )$values)

nn = 0:N
En = nn + 1/2 + 3*beta/8 * (nn^2 + nn + 1/4)
En


# alle EW's
plot(nn,res)
points(nn,En,col="red")


# first k EW's
k = 200
info = "exact En's (black)  vs.  Bohr-Sommerfeld En's (red),  beta = 10^(-2)"
plot(nn[1:k],En[1:k],col="red",main=info,xlab="n",ylab="En")
points(nn[1:k],res[1:k])
points(nn[1:k],En[1:k],col="red")