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?