# Non-numerical value for derivative at boundary using ParametricNDSolve to find eigenvalues using a simple shooting method

I have this code:

V[x_] = If[x <= -1, ∞, 1]*If[-1 < x < -a, 0, 1]*If[-a <= x <= a, V1, 1]*If[a < x < 1, 0, 1]*If[1 <= x, ∞, 1] /. {a -> 1/4}
sg = -psi''[x]/2 + V[x]*psi[x] == EE *psi[x] ;
sol = ParametricNDSolve[{sg, psi[-1]==0, psi'[-1]==1},psi, {x, -1, 1}, {EE, V1}]
Evaluate[psi[5, 30][1] /. sol]
Plot[Evaluate[psi[5, 30][x] /. sol], {x, -1, 1}]


I can't not explain the error:

ParametricNDSolve::ndnum: Encountered non-numerical value for a derivative at x == -1.. >>

and

Infinity::indet: Indeterminate expression 0. \[Infinity] encountered. >>

• I don't think you can use Infinity in your potential function V. You should instead specify the boundary conditions at x = -1 and x == 1 that the function should be zero. Commented Dec 8, 2015 at 21:18
• I tried to change V and set the boundary conditions: sol = ParametricNDSolve[{sg, psi[-1] == 0, psi[1] == 0, psi'[-1] == 1}, psi, {x, -1, 1}, {EE, V1}] But this can't be solved: "NDSolve is not currently able to solve boundary value problems with \ discrete variables. " Commented Dec 8, 2015 at 21:32
• Thank you for the accept! Usually it's a good idea to wait a day or so before accepting an answer for the purposes of attracting other answers (an accepted answer gets fewer looks). When you have enough reputation, please upvote the answers that you get that you like! Welcome to mathematica.SE. (I changed the title a bit; if you really feel like it, feel free to change it back, but it's a good idea for the titles to be as specific as possible). Commented Dec 9, 2015 at 0:14

## Original Post

The problem is that you have introduced a discontinuity in the derivative of your function by specifying that the potential energy for x <= -1 is infinite. If you remove the factors of If[x <= -1, ∞, 1] and If[x >= 1, ∞, 1], your code runs fine, because it now has no problem applying the boundary conditions at x == -1.

Here is your code, a little cleaned up. You can use one or the other of the definitions of v below:

v[x_, a_?Positive] := Piecewise[{{v1, -a <= x <= a}}]
v[x_, a_?Positive] := v1 UnitBox[x/(2 a)]
Plot[v[x, 1/4] /. v1 -> 10, {x, -1, 1}, Exclusions -> False]


Then, the Schrodinger equation is given by

sg[a_?Positive] := -psi''[x]/2 + v[x, a]*psi[x] == EE*psi[x];


and we solve this via

sol = ParametricNDSolve[{sg[1/4], psi[-1] == 0, psi'[-1] == 1}, psi, {x, -1, 1}, {EE, v1}];


Finally, we use FindRoot to see where the function goes to zero one left boundary and plot. For instance,

pot = 10;
energy = FindRoot[Evaluate[psi[e, pot][1] == 0 /. sol], {e, 5}]
(* { e -> 5.64844} *)


and

Plot[Evaluate[psi[e /. energy, pot][x] /. sol], {x, -1, 1}]


## Update

Let's see how well ParametricNDSolve does with this problem. The energy eigenvalues for the v1 -> 1 potential (i.e. just an infinite square well of width 2 units) should be

((π^2)/8)*Range[9]^2


My favorite way of quickly extracting numerical roots is as follows. Plot the function whose roots give us the energies with the roots found, i.e.

p1 = Plot[psi[e, 0][1] /. sol
, {e, 0, 100}, PlotRange -> {-0.5, 1}
, MeshFunctions -> Function[#2], Mesh -> {{0}}, MeshStyle -> {PointSize[0.015], Red}]


Extract the zeroes and refine them with FindRoot (this second step isn't entirely necessary, but it gets us better precision) via

rawEnergies = Sort@Cases[Normal@p1, Point[{a_, _}] :> a, Infinity]
refinedEnergies = e /. FindRoot[Evaluate[psi[e, 0][1] == 0 /. sol], {e, #}] & /@ rawEnergies


Finally, how close are they to the real values? Multiply by 8/π^2:

8/π^2 # & /@ refinedEnergies
(* {1., 4., 9., 16., 25., 36., 49., 64., 81.} *)


Excellent!

Finally:

Plot[Evaluate[psi[#, 0][x] & /@ refinedEnergies /. sol], {x, -1, 1}]


Those look like very nice sinusoidal functions!

Energies as a function of the height of the bump

Let's automate the procedure!

energies[v1_, maxE_] := Module[{p1, rawEnergies, refinedEnergies}
, p1 = Plot[psi[e, v1][1] /. sol, {e, 0, maxE}, PlotRange -> {-0.5, 1}, MeshFunctions -> Function[#2], Mesh -> {{0}}]
; rawEnergies = Sort@Cases[Normal@p1, Point[{a_, _}] :> a, Infinity]
; refinedEnergies = {v1, e /. FindRoot[Evaluate[psi[e, v1][1] == 0 /. sol], {e, #}]} & /@rawEnergies
]


This is somewhat inefficient, but it works well enough. We do this root-finding for a number v1s via

egys = Table[energies[v1, 100], {v1, 0, 30, 0.5}];


This spits out some errors because sometimes FindRoot has some problems, but here it doesn't matter. Finally, we plot this using

Show[
ListLinePlot[Flatten[egys, {2}]
, PlotRange -> {0, 40}
, Frame -> True, FrameLabel -> {"V", Style["E", Italic]}
, BaseStyle -> {FontFamily -> "CMU Serif", FontSize -> 15}
]
, Plot[Evaluate[\[Pi]^2/8 Range[9]^2], {x, 0, 30}, PlotStyle -> Dashed]
]


resulting in

The dashed curves are the particle-in-a-box energies. (By the way, the list of energies that we make might be ragged, because the energies get pushed up by the potential bump. For this reason, we have to transpose a ragged array, and we do this using the magic that is the second argument to Flatten.)

• Thanks a lot! It works just fine. Commented Dec 8, 2015 at 23:20
• But wait, there's more! Commented Dec 8, 2015 at 23:21