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:

 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*ω;
      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?

  • 1
    $\begingroup$ You should include your constants' values (eV, nm, etc) $\endgroup$ – Dr. belisarius Sep 25 '14 at 7:27
  • $\begingroup$ 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 $\endgroup$ – Ajay Sep 25 '14 at 7:34
  • $\begingroup$ 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 $\endgroup$ – Dr. belisarius Sep 25 '14 at 8:00
  • $\begingroup$ 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. $\endgroup$ – Ajay Sep 25 '14 at 8:09
  • $\begingroup$ Glad to help. Good luck with your assignment! $\endgroup$ – 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}]

Mathematica graphics

  • $\begingroup$ BTW I learned about it here mathematica.stackexchange.com/a/56164/193 $\endgroup$ – Dr. belisarius Sep 25 '14 at 7:55
  • $\begingroup$ 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) $\endgroup$ – Michael E2 Sep 25 '14 at 12:42
  • $\begingroup$ I meant learnt. Sometimes my English ... $\endgroup$ – Dr. belisarius Sep 25 '14 at 14:42

Not the answer you're looking for? Browse other questions tagged or ask your own question.