[?] Join solutions to the Scrödinger Equation for a barrier potential?

Questions
,
1

I want to obtain the the solutions to the time independent Schrödinger Equation for a barrier potential. I have created three equations:

tise1 = -\[HBar]^2/2 m D[\[Psi]1[x], {x, 2}] == e0 \[Psi]1[x]
tise2 = -\[HBar]^2/2 m D[\[Psi]2[x], {x, 2}] == (e0 - v0) \[Psi]2[x]
tise3 = -\[HBar]^2/2 m D[\[Psi]3[x], {x, 2}] == e0 \[Psi]3[x]

and the boundary conditions

bc = {\[Psi]1[0] == \[Psi]2[0], \[Psi]2[a] == \[Psi]3[a], 
D[\[Psi]1[x], x] == D[\[Psi]2[x], x] /. x -> 0, 
D[\[Psi]2[x], x] == D[\[Psi]3[x], x] /. x -> a}; ic = {}

(I intended to use ic to remove “left moving” waves right of the potential (x>a), but I could no figure out how to do this)
Next I used

sol = DSolve[
  Join[{tise1, tise2, tise3}, bc, ic], {\[Psi]1, \[Psi]2, \[Psi]3}, 
  x]

But the results are overwhelming and as far as I can tell: wrong. I would like to be able to get the wave function as an piecewise concatenation of psi 1-3 to add the time development and to calculate transmission and reflection rates.

Hopefully this makes sense. Any help (or pointers to such) will be appreciated.

/Mogens

2

First of all I think you should write

tise1 = -\[HBar]^2/(2 m ) D[\[Psi]1[x], {x, 2}] == e0 \[Psi]1[x]
tise2 = -\[HBar]^2/(2 m ) D[\[Psi]2[x], {x, 2}] == (e0 - v0) \[Psi]2[x]
tise3 = -\[HBar]^2/(2 m) D[\[Psi]3[x], {x, 2}] == e0 \[Psi]3[x]

Then, if e0 and v0 are contstants the solutions of the deqs are linear combinations Sin and Cos - functions yielding 6 constants of integration. You should try to fix these to get a solution to your problem.

3

I think it would be better to use DSolve on each piece separately and then apply the boundary conditions, rather than building them into one DSolve call.

4

This is just general advice. It is a good idea to rescale your problem, so there is no \[HBar] etc. at all. Also do not use Greek letters, code looks quite unreadable with them.

5

Frank,

Thanks for the advice, but when I solve them individually I get the same constants of integration for each solution (C[1] and C[2]). I there a way to avoid that, or do I need to replace them “by hand”. I suppose something like /. C[1] → c1 etc.

/Mogens

6

The documentation for DSolve has this example

In[1]:= DSolve[y''[x] == y[x], y[x], x, GeneratedParameters -> d]

Out[1]= {{y[x] -> E^x d[1] + E^-x d[2]}}
7

Thanks again, that is most helpful.

/Mogens