programming - 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?