# NDSolve, Schrödinger equation, and decaying solution

I am trying to solve a Schrödinger equation for a particle hitting a step potential using NDSolve in Mathematica. Here is my code:

mu = 6.;
m = mu;
R  = 5.;
Vs2 = 4./(2*m*R^2);
Vs = -10./(2*m*R^2) + Vs2;
Energy = 0.001
VCC[r_] = Vs*UnitStep[R - r] + Vs2*UnitStep[r - R];
L = 0;

system = {RC''[r] +
2/r*RC'[r] + (-L*(L + 1)/r^2 - 2*mu*(VCC[r] - Energy))*RC[r] ==
0, RC[0.001] == 1.0, RC'[0.001] == 0.0 };

syssol = NDSolve[system, {RC[r]}, { r, 0.001, 1000.}, MaxSteps -> 10000000];

Plot[Evaluate[{RC[r]} /. syssol], {r, 0.001, 200.0}, PlotRange -> {-1.1, 1.1}]


There should be decaying wave when particle hits the potential step, but NDSolve gives an increasing result. I am sure there is some trick to fix this, so I am waiting for you help.

• Before posting this question, you should first have followed the advice in the comments to your identical question on StackOverflow: stackoverflow.com/questions/9855619/… – Jens Mar 25 '12 at 18:38
• What is your m ? Also take a look at this (it's a bit different case but you can see how everything is set up): demonstrations.wolfram.com/ScatteringOverPotentialStep – Vitaliy Kaurov Mar 25 '12 at 18:39
• By the way, this has exact solutions, so I'm guessing this could be a homework problem and hence not appropriate for this forum. – Jens Mar 25 '12 at 18:43
• you forgot to define m here (this is to make it simpler for someone to cut and paste this code directly) – acl Mar 25 '12 at 18:43
• and, the root of you problem is mathematical, not an error in programming. did you try to do it analytically, for instance? – acl Mar 25 '12 at 18:44

noeckel’s answer on StackOverflow is spot on. This is not a Mathematica issue, this is a mathematical issue. Namely, Mathematica is giving you the correct solution to the system of differential equation and boundary conditions given. The conditions given (and in particular the derivative imposed at the origin) are incompatible with the expected decay. Bear in mind that, at $r \geq 5$, your wavefunction will have two components of the form $\exp(\alpha r)$ and $\exp(-\alpha r)$. For each set of boundary conditions, you get a different linear combination of these two, and the only conditions that make sense are those for which the diverging term is zero.
• +1. If memory serves me well, usual conditions for the radial wave function are the absence of increasing exponent at infinity and finiteness at zero (or at least, convergence of the integral of Abs[f[r]]^2 over the volume around zero. The divergence of this integral would correspond to the classical fall-on-the-center scenario, for potentials ~ 1/r^a with a > 2 IIRC). This rules out fixing the derivative at zero, as you said, and fixing the w.f. value itself does not make sense either, since it is fixed by w.f. normalization (that is, unless I forgot everything, it's been a while :)) – Leonid Shifrin Mar 25 '12 at 19:24