Schrodinger Equation by Shooting Method

Questions
,
1

Hello everybody

I want to solve the radial part of Schrodinger equation by shooting method. In this method I set the initial value of energy and appropriate boundary conditions. By varying the value of energy the best wave function, which satisfies the boundary conditions, is found. My code is:

ClearAll;
ClearSystemCache;
v[x_]:=If[0<x< 10,0,100];
energy=13.3;
xMax=12;
mc2=0.511 * 10^9;
hc = 1240 * 10^3 ;
solution = NDSolve[{x psi''[x]+ 2 psi'[x]== (8 Pi2 (mc2))/hc2 x ( v[x] - energy) psi[x], psi[0.01] == 1, psi'[0.01] == 0.001}, psi, {x, 0.01, xMax}];
Plot[psi[x]/. solution, {x, 0.01, xMax}]

Now I want to make a loop to change the value of energy and again recalculate the wave function until the correct wave function is obtained. I used Do[expr,{i,imin,imax,di}] but it didn’t work. Please guide me how to make a loop to change the energy and plot wave functions for each value of energy.

2

I think it would be better to use FindRoot .

3

Here’s an example of how it could be done. I wrote this before ParametricNDSolve was available.

In[1]:= shootsoln[d_?NumberQ] := 
 NDSolve[{y''[x] - 5 Sinh[5 y[x]] == 0, y[0] == 0, y'[0] == d}, 
   y, {x, 0, 1}][[1]]

In[2]:= yd[d_?NumberQ] := y[1] /. shootsoln[d]

In[4]:= soln = y /. (shootsoln[d] /.
    		FindRoot[yd[d] == 1, {d, -.05, .05}][[1]])

Out[4]= InterpolatingFunction[{{0., 1.}}, <>]
4

Dear Dr. Kampas
Thanks for your reply. My new code is:

CearAll;
ClearSystemCache;
v[x_] := If[0 < x < 10, 0, 100];
xMax = 12;
mc2 = 0.511 * 10^9;
hc = 1240 * 10^3 ;
system [energy_] := {
psi’‘[x] + 2 psi’[x] == (8 Pi^2 (mc2))/
hc^2 x ( v - energy) psi, psi[0] == 1, psi’[0] == 0.001};
solution[energy_] = NDSolve[system [energy], psi[x], {x, 0, xMax}];
yend[energy_?NumericQD] := First[psi[x] /. solution [energy] /. x];
bc = FindRoot[yend [energy] == 0, {energy, 13, 13.5}];
Plot[ Evaluate[psi[x] /. solution[energy /. bc]], {x, 0, xMax}]

?again it doesn’t give the correct answer. it should be something like the attached file. one more thing is that this approach doesn’t give us the energy value for which the wave function has been plotted.
What should I do?

5

It would be easier to read your code if you put it into a “code sample” block, which is the first icon.

I suggest that you first test your solution[energy] function to see if it gives the correct result for the correct energy. Also note that the code I posted as an illustration can be split up to get the value found by FindRoot. Since there are multiple eigenvalues, you need to give FindRoot starting values near the correct result.

In[1]:= shootsoln[d_?NumberQ] := 
 NDSolve[{y''[x] - 5 Sinh[5 y[x]] == 0, y[0] == 0, y'[0] == d}, 
   y, {x, 0, 1}][[1]]

In[2]:= yd[d_?NumberQ] := y[1] /. shootsoln[d]

In[3]:= fr = FindRoot[yd[d] == 1, {d, -.05, .05}]

Out[3]= {d -> 0.0457504}

In[4]:= soln = y /. (shootsoln[d] /. fr)

Out[4]= InterpolatingFunction[{{0., 1.}}, <>]
6

Dear Dr. Kampas

Thanks a million for all your hints and suggestions. I used ParametricNDSolve function and fortunately it worked. now I want to have legend on the graph to label wave function corresponding to each energy value. how should I do that? is that PlotLegend function?

7

The PlotLegends option for your plot function may do what you want.