Increasing NDSolve computation range for 3D stationary PDE

Question

I am trying to solve 3D Schrodinger equation (with mass equals to 1) and with no potential (I will add a potential later). My goal is to set boundary conditions on the outer boundary of a sphere and see that the result obtained by NDSolve is consistent with it.

The analytical solution for this equation is simply the spherical harmonics times the spherical Bessel function.

So first, I defined:

rCoor[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2]
θCoor[x_, y_, z_] := ArcTan[z, Sqrt[x^2 + y^2]]    
ϕCoor[x_, y_, z_] := ArcTan[x, y]
YlmCart[l_, m_, x_, y_, z_] := SphericalHarmonicY[l, m, θCoor[x, y, z], ϕCoor[x, y, z]]

Now I'm trying to solve the equation, with l=m=0 (the spherical harmonic paramters) and for energy parameter of Ef = 1:

Needs["NDSolve`FEM`"];
l = 0;
m = 0;
Ef = 1;
kf = Sqrt[2 Ef];
rmax = 5 kf^-1;
mesh = ToElementMesh[Ball[{0, 0, 0}, rmax]];

boundary[l_, m_, k_, x_, y_, z_] := 
 SphericalBesselJ[l, k rCoor[x, y, z]] YlmCart[l, m, x, y, z]
op = -(1/2) Laplacian[u[x, y, z], {x, y, z}] - Ef u[x, y, z];
sol = NDSolveValue[{op == 0,DirichletCondition[u[x, y, z] == 
    boundary[l, m, kf, x, y, z] , True]}, u, {x, y, z} ∈ mesh];

Now I'm plotting the result to compare with the analytical solution (checking only in the x direction, in which theta and phi are constant):

Plot[{Re[ sol[x, 0, 0]], 
  Re[boundary[l, m, kf, x, 0, 0]]}, {x, 0, rmax}]

as you can see, the result is just fine and perfectly aligned with the analytical result. Now, when I'm increasing the range of computation by changing rmax, the result start to deviate. This is the result for rmax = 30 kf^-1:

What is going on here? How can I solve this problem?

Answer 1

For rmax = 5 kf^-1, ToElementMesh discretizes the region of integration by default into about 7000 elements. For rmax = 30 kf^-1, it does the same. Thus, resolution is much reduced, and this is the source of the inaccuracy for the larger rmax. To increase resolution and display the mesh, use

mesh = ToElementMesh[Ball[{0, 0, 0}, rmax], MaxCellMeasure -> {"Length" -> 1.2}]
mesh["Wireframe"]

instead of the definition of mesh in the question. (It also is necessary to add Needs["NDSolveFEM"], which apparently was omitted from the question by oversight.) It yields a mesh of about 90,000 elements

which is as many as my computer can accommodate without becoming extremely slow. The solution then is

which is less accurate than desired but much better than the result depicted in the question. To obtain still better results, one should use some of the many options in ToElementMesh to generate a mesh with high resolution radially and low resolution azimuthally.