Harmonic oscillator solver

In this notebook we solve:

	-  (((d^2)ψ)/(dx^2)) + (1/2)(Subscript[V, 0])(x^2)= E ψ
	
in the region  -xmax < x < xmax.  We assume zero value boundary conditions for this problem.

xmax=10;  (* The system goes from -xmax to xmax*)
ndiv=10;  (* Number of intervals in each unit spacing *)
dx=N[1/ndiv] ; (* Our Δx *)
nmax=2*xmax*ndiv+1;  (* Total number of points solved *)
invdxsq=1/(dx^2)  ;   (* A useful constant *)

V0=1.;    (*  Strength of potential  *) 
V[x_]:=(V0/2)* (((x))^2);

Here’s a plot of the potential:

vdata=Table[x=(i-1)*dx-xmax;{x,V[x]},{i,1,nmax}];
ListPlot[vdata]

Clear[hmat];
hmat=Table[0,{i,1,nmax},{j,1,nmax}];
Do[hmat[[i,i]]= invdxsq +V[(i-1)*dx-xmax],{i,1,nmax}];
Do[hmat[[i,i+1]]=-invdxsq /2,{i,1,nmax-1}];
Do[hmat[[i+1,i]]=-invdxsq /2,{i,1,nmax-1}];

eigenout=Eigensystem[hmat];
eigval=eigenout[[1]];
eigveclist=eigenout[[2]];
sortedeigs=Sort[eigval];

The eigenvalues are not sorted.  We can sort them if we wish:

ListPlot[sortedeigs]

Note that they are linear for small energies.  They should be equal to  (n+(1/2)).  They are pretty close for the first few.

Take[sortedeigs,8]

This function picks out the i-th eigenenergy in our original list.

findIndex[i_]:=Position[eigval,sortedeigs[[i]]]

Let’s find the index (location) of the lowest eigenenergy

findIndex[1]

Ok, now we plot tht 201st eigenvector:

ListPlot[eigveclist[[201]],PlotRange-> All]

And here is the second :

findIndex[2]

ListPlot[eigveclist[[200]],PlotRange-> All]