8

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?

xzczd
  • 65,995
  • 9
  • 163
  • 468
Ajay
  • 227
  • 2
  • 9
  • 1
    You should include your constants' values (eV, nm, etc) – Dr. belisarius Sep 25 '14 at 07: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 07: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 08: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 08:09
  • Glad to help. Good luck with your assignment! – Dr. belisarius Sep 25 '14 at 08:16

1 Answers1

9
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

Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453