ParametricNDSolveValue for ODE with boundary at infinity

Question

I'm trying solve the Schrödinger equation for a given potential using the function ParametricNDSolveValue following the the first method of the post:

Find eigen energies of time-independent Schrödinger equation

It works fine, but it bothers me that we have to impose the somewhat artificial boundary condition ψ[1]==1 instead of the real condition ψ[Infinity]==0.

I've tried two different ideas to get around this:

a. Imposing the condition ψ[100]==10^-5 or similar values, the problem with this is that afterwards Findroot doesn't return the right values of the energy (I tested with the known results. of the harmonic oscillator). Oddly, if I use ψ[3]==0.01 I get better results than with ψ[3]==0.0001 or ψ[10]==0.0001 and definitely than with ψ[100]==10^-5 (?!)

b. Using ParametricNDSolveValue with two parameters, one for the Energy and one for the the boundary conditions, ie. ψ[100]==k, and the using Findroot for the energy making k=0, but again, results are wrong.

So my questions are:

1. Why do I get better results using the boundary condition ψ[3]==0.01 than using ψ[3]==0.0001 or ψ[10]==0.0001 when the latest are closer to the real boundary condition?

2. Is there a way of using the fonction ParametricNDSolveValue when boundaries are infinity or close enough?

Answer 1

1. Why do I get better results using the boundary condition ψ[3]==0.01

E^(-(x)^2/2)*(1/Pi)^(1/4) /. x -> 3 // N gives 0.00834425, so to answer your question: ψ[3]==0.01 gives a better result because it is a better approximation, for the fundamental state at least. And ψ[3]==0.008 would be even better!

But you should be aware that this will not be a good approximation for higher energy states.

2. Is there a way of using the fonction ParametricNDSolveValue when boundaries are infinity or close enough?

Since computers dont have infinite memory we cant solve DE over infinitely large regions. We have to truncate our region of intrest at some point.

I found that

U[x_] := 1/2  x^2;
L = 6;
pfun = ParametricNDSolveValue[{-1/2 \[Psi]''[x] + U[x] \[Psi][x] == 
     Ei \[Psi][x], \[Psi][L] == 0, \[Psi][0] == 
     1}, \[Psi], {x, -L, L}, {Ei}];

Plot[Evaluate@pfun[.5][x], {x, -3, 3}, 
 PlotRange -> {All, {-2, 2}}]

Plot[Evaluate@pfun[2.5][x], {x, -3, 3}, 
 PlotRange -> {All, {-2, 2}}]

works well for even half energies (I tried up to 6.5). If you want more excited states you might need to truncate over a larger area (bigger L), but that will make the computation take more time.

Also notice the initial condition \[Psi][0] == 1, it was necessary to keep the final result from being the trivial one: \[Psi][x] == 0, for odd solutions you will have to change it to \[Psi]'[0] == 1. Both artificial conditions might generate unnormalized solutions.