# Eigen energies of 2D time independent Schrödinger equation

Firstly, I am aware that a similar question of this type was asked, but it was not helpful as it was only one dimensional (so please don't mark my question as a duplicate)

The trouble I am having is using ParametricNDSolveValue to solve my Schrödinger PDE. Here is what I have:

ℏ = 1;
m = 1;
Δ = 10^-3;
k = 1;
X = 70;
V[x_, y_] := -k/Sqrt[1 + x^2 + y^2]
Sol = ParametricNDSolveValue[{-ℏ^2/(2 m) (D[ψ[x, y], {x, 2}] + D[ψ[x, y],
{y, 2}]) + V[x, y] ψ[x, y] == Etr ψ[x, y], ψ[x, -X] == Δ, ψ[x, X] == Δ,
ψ[-X, y] == Δ, ψ[X, y] == Δ}, ψ[0, 0], {x, -X, X}, {y, -X, X}, {Etr}]


It seems to work okay initially, but when I try to find roots of this function ie:

FindRoot[Sol, {-0.45}]


I get values that are inconsistent with my solution in polar/cylindrical coordinates. That solution gave me proper values and good solutions (because it was done in essentially 1 dimension) If you do the following plot:

EList = Range[-0.46, 0, 0.0005];
ETable = ParallelTable[{EList[[i]], Sol[EList[[i]]]}, {i, 1,
Length[EList]}];
ListPlot[ETable]


The result is terrible. Does anyone have any idea as to how to solve this issue, or whether it is a bug or not?

• By imposing fixed values of the function on a square you're forcing an unphysical solution to this rotationally symmetric problem. It will probably violate the decaying boundary condition at infinity. So the setup of the problem in this 2D calculation is most likely incorrect, but definitely inconsistent with any kind of solution in polar coordinates that you mention as a comparison case. – Jens Apr 17 '15 at 4:21

The calculation seems to work, producing the curve,

Plot[Sol[x], {x, -.45, 0}]


Finding zeroes of Sol requires a modification of the expression in the Question,

FindRoot[Sol[x], {x, -0.45}, Evaluated -> False]
(* {x -> -0.325184} *)


It is not possible to comment on why this might not agree with another calculation without seeing that other calculation.

Your code has two problems: ψ[0, 0] should be ψ inside your ParametricNDSolveValue, and your boundary conditions are not physical.

Since I dont have your analitical solutions here I tried to guess some good boundary conditions and came up with Exp[-Sqrt@(r^2)], which has the necessary symmetry.

Also, I reduced the interval of integration since 70 seemed excessive. If I am wrong it is easy to fix. And to make the energy have positive values for bound states I lifted your potential by adding k to it.

The code is as follows:

ClearAll["Global*"];
\[HBar] = 1;
m = 1;
\[CapitalDelta] = 10^-3;
k = 1;
X = 5;
V[x_, y_] := k (1 - 1/Sqrt[1 + x^2 + y^2]);
ini[x_, y_] := Exp[-Sqrt@(x^2 + y^2)];
bc = \[Psi][x, y] ==
ini[x, y] /. {{x -> X}, {x -> -X}, {y -> X}, {y -> -X}};
sys = {-\[HBar]^2/(2 m) (D[\[Psi][x, y], {x, 2}] +
D[\[Psi][x, y], {y, 2}]) + V[x, y] \[Psi][x, y] ==
Etr \[Psi][x, y], bc};
sol = ParametricNDSolveValue[
sys, \[Psi], {x, -X, X}, {y, -X, X}, {Etr}]


now to find the eigen energy of the fundamental state by trial and error, I found that Etr==0.569717 gives a nice looking function. The closer you are to the real value the better behaved the solution is. So the absence of divergences is a good sign.

Plot3D[sol[.569717][x, y], {x, -X, X}, {y, -X, X}, PlotRange -> All]
`