(Edited) For finding the ground state wave function of:

$ H\psi(x) = (-1/2)d^2\psi(x)/dx^2 + (1/2)x^2\psi(x) = E \psi(x)$

I have written:

mOneDSchEq[n_] :=
Table[Switch[i - j, -1, p[x[i]], 
0, (10/(n + 1))^2 q[x[i]] - 2 p[x[i]], 1, p[x[i]], _, 0], {i, 
n}, {j, n}];

q[x_] := -x^2; p[x_] := 1;
Xarray[n_] := Do[x[i] = -5 + i 10/(n + 1), {i, 0, n + 1}];

EigVec[n_] := Eigenvectors[mOneDSchEq[n]];
lisEigVec = EigVec[35];
OneEigVec[j_] := Part[Reverse[lisEigVec], j];
y[i_] := Part[OneEigVec[1], i];
listOfPoints = 
Join[{{x[0], 0}}, Table[{x[i], y[i]}, {i, 1, 35}], {{x[36], 0}}];
ListPlot[listOfPoints, PlotJoined -> True, PlotRange -> All, 
PlotLabel -> "Ground State Wave Function of Harmonic Oscillator", 
AxesLabel -> {"x", "y"}]

Which I have obtained the Gaussian, correctly.

The question that came to my mind is that:

Is it possible by knowing the ground stat eigenvalue, i.e. 1/2, solving the Schrödinger equation numerically, and obtain the ground state wave function? in other words, to solve:

$ H\psi(x) = (-1/2)d^2\psi(x)/dx^2 + (1/2)x^2\psi(x) = (1/2) \psi(x)$


$ H\psi(x) = (-1/2)d^2\psi(x)/dx^2 + (1/2)x^2\psi(x) = (3/2) \psi(x)$

So, I wrote:

s = NDSolve[{-(1/2) \[Psi]''[x] + (1/2) x^2(\[Psi][x]) == (1/2) \[Psi][
  x], \[Psi][-5] == 0, \[Psi][5] == 0}, \[Psi], {x, -5, 5}]

Plot[Evaluate[\[Psi][x] /. s], {x, -5, 5}, PlotRange -> All]

BUT, I got nothing. What is the problem?

The other question is that, I was traveling through the website and found an elegant approach to two dimensional Harmonic Oscillator here.

My question is, if again we want to solve the Schrödinger equation numercally and obtain wave functions, now two dimensional, by knowing the eigenvalues, what should we do? For example:

$ H\psi(x,y) = (-1/2)(d^2/x^2 + d^2/dy^2)\psi(x,y) + (1/2)(x^2 + y^2)\psi(x,y) = (1) \psi(x,y) $


$ H\psi(x,y) = (-1/2)(d^2/x^2 + d^2/dy^2)\psi(x,y) + (1/2)(x^2 + y^2+ x y)\psi(x,y) = (0.96) \psi(x,y) $

Thanks for your attention!

  • $\begingroup$ In Jens' answer, isn't the 1/(2 a^2) bit there to take into account the factor of 1/2 in front of the laplacian? Also, the Partition is there because he is representing 2d space in a 1d vector (basically, he discretises space, then take the 2d matrix and set the rows one after the other to each other so as to form a 1d vector; the Partition undoes this). I'll write an answer tomorrow if I remember, have time and nobody else does! $\endgroup$
    – acl
    Jun 23, 2014 at 0:16
  • $\begingroup$ Actually I removed "2" in the denominator of 1/(2 a^2), but the answer seemed to me non-sense! and also, why did not plot the two dimensional wave function directly, and he used ListDensityPlot? $\endgroup$
    – SomeBody
    Jun 23, 2014 at 0:21
  • $\begingroup$ Hmm sorry, I didn't understand; can you rephrase? $\endgroup$
    – acl
    Jun 23, 2014 at 0:24
  • $\begingroup$ I removed 2 in the denominator and also 1/2 in the potentials, so I expect to obtain two times of energy, i.e. 2E; but It was not. $\endgroup$
    – SomeBody
    Jun 23, 2014 at 0:26
  • $\begingroup$ to Jens: I tried to comment there, but I was not allowed! $\endgroup$
    – SomeBody
    Jun 23, 2014 at 0:28

1 Answer 1


To give another answer for the one-dimensional harmonic oscillator, let's use a different approach based on the NDSolve functionality I alluded to in the linked answer. Edit: I also update the linked answer to include the analogue of this approach in two dimensions.

n = 2000;
a = .02;
grid = N[a Range[-n, n]];
derivative2 = 
 NDSolve`FiniteDifferenceDerivative[2, grid]["DifferentiationMatrix"]


potential = Map[(1/2 #^2) &, grid];

hamiltonian = -derivative2/2 + 

eigenvalues = Chop[Eigenvalues[hamiltonian, -10]]

{9.5, 8.5, 7.5, 6.5, 5.5, 4.5, 3.5, 2.5, 1.5, 0.5}

v = Chop[Eigenvectors[hamiltonian, -10]];

ListLinePlot[{Abs[v[[-1]]]^2, Abs[v[[-2]]]^2, 
  Abs[v[[-3]]]^2}, DataRange -> grid[[{1, -1}]], 
 PlotRange -> {{-4, 4}, All}]

plot of three functions

Here I used a grid spacing of a = 0.02 and get numerically very exact solutions for the lowest states of the harmonic oscillator.

The matrix representing the second derivatives (derivative2) in the Laplacian is generated using FiniteDifferenceDerivative.

To address some of the other issues in the question:

The initial code in the question didn't produce a result for me. However, since you state you got the desired result, I assume that there is some typo in the question. Definitely, one can improve the first code block by wrapping the generated Hamiltonian matrix in N to make it into a machine-precision matrix that can be diagonalized much faster.

However, the main question seems to have been: why does the differential equation

s = 
 NDSolve[{-(1/2) ψ''[x] + (1/2) x^2 (ψ[x]) == (1/2) ψ[
      x], ψ[-5] == 0, ψ[5] == 0}, ψ[x], {x, -5, 5}];
ψSol[x_] = ψ[x] /. s[[1, 1]];

Plot[Evaluate[ψSol[x]], {x, -5, 5}, 
 PlotRange -> All]

yield an apparently empty plot? The answer is that the boundary conditions are incorrect if you're looking for a non-trivial solution. The solver actually finds the only possible answer, $\psi(x)\equiv 0$ for all $x$. But this is because you forced the wave function to be zero at two points whereas the ground state by definition has no nodes!

So you should solve the following equation instead:

s = 
 NDSolve[{-(1/2) ψ''[x] + (1/2) x^2 (ψ[x]) == (1/2) ψ[
      x], ψ[0] == 1, ψ'[0] == 0}, ψ[x], {x, -5, 5}];
ψSol[x_] = ψ[x] /. s[[1, 1]];

Plot[Evaluate[ψSol[x]], {x, -5, 5}, 
 PlotRange -> All]


This yields the expected Gaussian. I chose boundary conditions for the function to be 1 and its derivative to be 0 at the origin.

  • $\begingroup$ It is much more complicated than I expected. But, as I asked in the question, if, for example the eigenvalue is given 1/2, why does not my simple approach in the question work? $\endgroup$
    – SomeBody
    Jun 23, 2014 at 1:21
  • $\begingroup$ @JensI assumed -5 and 5 physical infinity, where the wave function vanishes. $\endgroup$
    – SomeBody
    Jun 23, 2014 at 1:42
  • $\begingroup$ But then you can't blame Mathematica for giving you the exact solution satisfying your approximate assumption. $\endgroup$
    – Jens
    Jun 23, 2014 at 1:44
  • $\begingroup$ And in your two dimensional isotropic case, if E = 1 is given for example, is your plot there, reproducible by solving the differential equation for E = 1? $\endgroup$
    – SomeBody
    Jun 23, 2014 at 1:47
  • $\begingroup$ No, that won't work. The real problem is that this is an elliptical boundary-value problem, for which NDSolve is not always able to give the solutions you want. That's why the approach based on matrix diagonalization is still necessary. In one dimension, the boundary-value problem isn't too hard because there is only one independent variable. But when there are more than 2 variables, you have a partial differential equation for which NDSolve is usually not able to find a solution. $\endgroup$
    – Jens
    Jun 23, 2014 at 1:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.