Solving for the eigenvectors and eigenvalues of the quantum hamiltonian

quantum simulation.ma

(*Define the potential energy function*)
potential[x_] := x^2

(*Set the mass and Planck's constant*)
mass = 1;
hbar = 1;

(*Set the range of x values and the number of points to use*)
xMin = -5;
xMax = 5;
numPoints = 1000;

(*Calculate the step size*)
stepSize = (xMax - xMin)/numPoints;

(*Generate the list of x values*)
xValues = Range[xMin, xMax, stepSize];

(*Calculate the kinetic energy operator matrix*)
kineticMatrix = 
  Table[If[i == j, -hbar^2/(2*mass*stepSize^2), 
    If[Abs[i - j] == 1, hbar^2/(2*mass*stepSize^2), 0]], {i, 
    numPoints}, {j, numPoints}];

(*Calculate the potential energy operator matrix*)
potentialMatrix = Table[potential[xValues[[i]]], {i, numPoints}];

(*Calculate the Hamiltonian matrix*)
hamiltonianMatrix = kineticMatrix + DiagonalMatrix[potentialMatrix];

(*Solve for the eigenvalues and eigenvectors of the Hamiltonian \
matrix*)
{eigenvalues, eigenvectors} = Eigensystem[hamiltonianMatrix];

(*Select the ground state wave function*)
groundStateWaveFunction = eigenvectors[[All, 1]];

(*Normalize the wave function*)
groundStateWaveFunction = 
  groundStateWaveFunction/Norm[groundStateWaveFunction];

(*Plot the wave function*)
ListLinePlot[groundStateWaveFunction, DataRange -> {xMin, xMax}, 
 PlotRange -> All]