Calculating eigenvectors/values from Hamiltonian

[fenny]

Homework Statement

I've constructed a 3D grid of n points in each direction (x, y, z; cube) and calculated the potential at each point.
For reference, the potential roughly looks like the harmonic oscillator: V≈r²+V₀, referenced from the center of the cube.
I'm then constructing the Hamiltonian by using a 2^nd order numerical derivative stencil in each direction representing the Laplacian seen in the TISE, multiplying each stencil by the factor $\frac {ħ} {2m}$, and adding the potential.

Homework Equations

$H = \frac {-ħ^2} {2m} ∇^2 + V$

$∇^2 = \frac {∂^2} {∂x^2} +\frac {∂^2} {∂y^2} +\frac {∂^2} {∂z^2}$

The Attempt at a Solution

My approach was to simply solve for the eigenvalues and corresponding eigenvectors of H (currently using "eig" function in MATLAB). However, I generally observe peaks at the corners or edges when verifying the probability density observed ($Ψ* Ψ$) and approximately zero elsewhere in the cube.
I've gone through my equations and verified the units to be correct a few times now, and I'm wondering if something is being missed on my end or if something is fundamentally flawed with this approach. Eventually I would also like to solve the problem using the variational method as well to compare the answers (and double check my calculation/method).
For a grid with n points along each axis, n³ eigenvectors/values are obtained due to construction of a n³ × n³ matrix representing H. I initially examined the lowest energy states primarily but am now sampling along the different energies and still observe the same behavior.
I'm also wondering if my potential (V) should be referenced to the lowest point (near the center), essentially removing the constant term. My fear was that the increase along the H matrix diagonal could be skewing a calculation. As it is, my potential terms (added to the primary diagonal) appear to be ~10^-10 while my stencil terms considering the neighboring points (adjacent to Ψ_i,j,k) are approximately ~10^-15 due to small reduced Planck constant.

 

[hilbert2]

You need to do some kind of change of units that makes the terms closer to unity, for example $\hbar = 1, m=1$. Also, make sure that the matrix eigenvalue equation $H\psi = E\psi$ has the factors of $1/(\Delta x)^2$ in the matrix, not in the RHS vector (the nonzero matrix elements will then become relatively large numbers). That way you should get a decent solution.