# NDSolve and Piecewise functions [duplicate]

Backslide introduced in 9.0, persisting through 11.0.1

I want to clarify upfront, This is not a please help me get the answer to my homework question. I know what I need to do I just do not know how to get mathematica to do it.

I am working on an assignment from my Quantum Mechanics professor, where we attempt to glean the form of the wave functions for the first two eigenstates of a single electron in a potential. We are just taking a guess at the energy and using NDsolve to get a solution to the problem. Then inspecting it visually to make sure it looks right. Here is what I have for part I the harmonic Potential:

 (*Constants*)
me = 9.11*10^(-31);
h = 6.62607 10^(-34);
hbar = h/(2 Pi);
nm = 1.00*10^(-9);
eV = 1.60218 10^(-19);
U = 1.5*eV;
L = 1.5*nm;

k = eV/nm^2;
Vhar[x_] := (k*x^2)/2;
ω = (k/me)^(1/2);
L = 1.5*nm;
α = Sqrt[(me*k/hbar^2)];
Ehar = (ν + 1/2) hbar*ω;
Do[
Print[
Plot[Evaluate[
f[x] /. NDSolve[{(-hbar^2/(2 *me))
f''[x] == (((ν + 1/2) hbar*ω) - Vhar[x]) f[x],
f[L] == 0, f'[-L] == 0.1}, f, {x, -L, L}]], {x, -L, L},
PlotRange -> All, PlotStyle -> Thick]], {ν, 0, 1, 1}];


Works Flawlessly, Great! I have a similar solution for the particle in a rectangular potential with an infinite wall at X=0, and a Finite one at x=L.

The last part of the problem wants us to look at a potential with infinite walls at x=+/=L, and a barrier of height U at -0.1L

Here I cant get the equations to solve. I have the Potential for part two in the form of a piece wise function but there appears to be an issue with the new one for part 3. My code is below:

Vb3[x_] := If[x < -L, 10^10, If[-0.1*L < x < 0.1*L, U, 0]];
NDSolve[{(-hbar^2/(2 *me)) f''[x] == (2.40986*^-20 - Vb3[x]) f[x],
f[L] == 0, f'[-L] == 0.1}, f, {x, -L, L}];


I get the following error.

NDSolve is not currently able to solve boundary value problems with discrete variables.

I cant find anything in the docs about how to fix this, but I know the problem is able to be solved. So how can I change my setup to avoid this problem with the boundary condition error?

• You should include your constants' values (eV, nm, etc) – Dr. belisarius Sep 25 '14 at 7:27
• They are just unit conversions, and don't really apply to the question. I wanted to be as brief as possible while getting the questions across. They have been added – Ajay Sep 25 '14 at 7:34
• Searching for your error message in this site you can find many related problems, since square (or triangular or whatever) potentials are the cool thing in QM – Dr. belisarius Sep 25 '14 at 8:00
• It sounds silly now saying it out-loud but never tried that. Thank you for your help. That question you linked has a ton of information to help me make this assignment looks a lot nicer. – Ajay Sep 25 '14 at 8:09
• Glad to help. Good luck with your assignment! – Dr. belisarius Sep 25 '14 at 8:16

Vb3[x_] := If[x < -L, 10^10, If[-0.1*L < x < 0.1*L, U, 0]];

u = NDSolve[{(-hbar^2/(2*me)) f''[x] == (2.40986*^-20 - Vb3[x]) f[x],
f[L] == 0, f'[-L] == 0.1}, f, {x, -L, L},
Method -> {"DiscontinuityProcessing" -> False}]
Plot[f[k] /. u, {k, -L, L}] • BTW I learned about it here mathematica.stackexchange.com/a/56164/193 – Dr. belisarius Sep 25 '14 at 7:55
• Hmm, the other solution in that answer, essentially Method -> "FiniteElement", produces almost the same picture, but off by a factor of about 1.65*^37 -- which is huge. But this one seems correct. (+1) – Michael E2 Sep 25 '14 at 12:42
• I meant learnt. Sometimes my English ... – Dr. belisarius Sep 25 '14 at 14:42