The subject of potential barriers has to do with finding the reflection and transmission properties of particles through potential walls. The energy is given , it is not an unknown. The basic solutions are either exponentials with imaginary arguments like exp(+ikx) ,

exp(-ikx) , or real exponentials exp(-kx) , exp(+kx) .

Our first example is that of a step potential of height V . Two cases will be considered , That of E < V where particles come from the left to the right and the case when E >V .

Fig 4.1

The energy is given , this is not an eigenvalue problem. We want to find the solutions Ψ in the two regions and the transmission and /or reflection characteristics. Schrodinger equation for each region is

d² Ψ₁ /dx² = -(2mE/h'² ) Ψ₁ , x ≤ 0 (1)

and d² Ψ₂ /dx² = {2m(V-E)/h'² } Ψ₂ , x ≥ 0 . (2)

The region x>0 is classically forbidden beacause there E < V and therefore the kinetic energy T = E-V would be negative.

It is required that at the boundary (x=0)

Ψ₁ (0) = Ψ₂(0) , (3)

(dΨ₁ /dx)_x=0 = (dΨ₂ /dx)_x=0 . (4)

Define two wave numbers k₁ = (2mE/h'² )^1/2 and k₂ = {2m(V-E)/h'² }^1/2 .

The general solutions of (1) and (3) are

Ψ₁ (x) = A exp( +ik₁x ) + B exp( -ik₁x ) (5)

If we wish to write the time dependent function , it is

Ψ₁ (x,t) = A exp( +i(k₁x-ωt) ) + B exp( -i(k₁x +ωt ) ) . (6)

The term ~ A is a plane wave traveling to the right and the term ~ B is travelling to the left.

Given that the probability current density is J = (ih'/(2m)){ Ψ ∂ Ψ^* /∂x - Ψ^* ∂ Ψ /∂x } (see eq.59 Chap 3) we apply it to Aexp(ik₁ x).

One obtains J= (ih'/(2m)) { A exp( ik₁ x) ( -ik₁) A^* exp( -ik₁ x) - A^* exp( -ik₁ x) (ik₁ )Aexp(ik₁ x) }

= ( hk₁/m) A A^* = (p/m) A A^* = v A A^* ~ 1/time (7)

Eq. (7) represents the flux of particles travelling to the right. So J_right = v A A^* and the flux of particles reflected and travelling to the left is obtained changing p→ -p in (7) , therefore

J_left = -v B B^* . (8)

For x > 0 the mathematical solutions are

Ψ₂ (x) = C exp( - k₂x ) + D exp( + k₂x ) . (9 )

One must set D = 0 otherwise the solution diverges as x → ∞ , leaving

Ψ₂ (x) = C exp( - k₂x ) . ( 10 )

The flux of Ψ₂ (x) is zero as seen immediately from

J = (ih'/(2m)){ Ψ ∂ Ψ^* /∂x - Ψ^* ∂ Ψ /∂x}

= (ih'/(2m)){ C exp( - k₂x )(-k ₂C^*exp exp(-k₂ x) ) - C^* exp(-k₂ x) (- k ₂Cexp exp(-k₂ x) ) }

= (ih'/(2m)) { - k₂ C C^* exp(-2k₂ x) + k₂ C C^* exp(-2k₂ x) } = 0 .

This means that while there is a probability of finding the particle in the region x>0 ( classically frobidden) , there is no net flux.

Another implication is that J_right + J_left = 0 that is, A A^* = B B^* .

From letting C= 1 and applying the conditions of continuity , eqs 3-4, it follows that

A + B = 1 , (11)

ik₁ (A - B) = -k₂ , (12)

whence

A = (1/2) ( 1 +i k₂ /k1) (13)

B= (1/2) (1 - ik₂ /k₁ ) . (14)

Notice that A= B* when E<V . Substitution of (13) , (14) in (5) gives ,

Ψ₁ (x) = cos(k₁x) - (k₂ /k₁) sin (k₁x) for x < 0 (15)

and Ψ₂ (x)= exp(-k₂ x) , x > 0 . (16)

Example 4.1 : Let V=1 , E= (1/2) V = 1/2 . Then k₁ = (2E)^1/2 = 1 , k₂ =(2(V-E))^1/2 = 1.

Taking another look at Ψ₁ (x) = cos(k₁x) - (k₂ /k₁) sin (k₁x) (a real function) we conclude that the current density is J=0.

Fig. 4-2 . Psi in the case in which E=V/2.

FORTRAN code for Fig 4.1

c plot of psi for potential barrier E < V
real k1,k2
data k1,k2/1.,1./
data xi,xf,nx/-7.,3.,70/
dx=(xf-xi)/float(nx)
do 10 i=0,nx
x=xi+dx*float(i)
print 100,x,psi(k1,k2,x)
10 continue
100 format('x,psi=',2(3x,e10.3))
stop
end

function psi(k1,k2,x)
real k1,k2
if(x.lt.0.)psi=cos(k1*x)-(k2 /k1)*sin(k1*x)
if(x.ge.0.)psi=exp(-k2*x)
return
end

Numerical solution of equations (1) and (2). The FORTRAN code given below integrates the DE using finite differences.

It starts at x= 0 with the initial values, psi0=1. and psi1=psi0 + dx*(-k2*psi0).

Then

psi2=2.*psi1-psi0 -dx**2*(2.*E*psi1) x < 0

and

psi2=2.*psi1-psi0 -dx**2*{2.*(E-V)*psi1} x>0 .

dx or ∆x has to be small relative to the Length scale which happens to be L_scale =(h'² /(mE) )^1/2 = 1/E^1/2 .

Fig 4-3. Numerical integration of Schrodinger's equation.

Fig 4-3 b . Plot of Psi² .

FORTRAN code for the numerical integration

c plot of psi for potential barrier E < V
real k1,k2
data E,V/.5,1./
data xi,xf,nx/0.,3.,1000/
k1=sqrt(2.*E)
k2=sqrt(2.*(V-E))
dx=(xf-xi)/float(nx)
kp=int(float(nx)/40)
kount=kp
psi0=1.
psi1=psi0+dx*(-k2*psi0)
print 100,xi,psi0
do 10 i=0,nx
x=xi+dx*float(i)
psi2=2.*psi1-psi0 -dx**2*d2(psi1,x-dx,E,V)
if(i.eq.kount)then
print 100,x,psi2
kount=kount+kp
endif
psi0=psi1
psi1=psi2
10 continue
100 format('x,psi=',2(3x,e10.3))
stop
end

function d2(psi,x,E,V)
if(x.lt.0.)d2= 2.*E*psi
if(x.ge.0.)d2=2.*(E-V)*psi
return
end

The coefficients A , B can be determined from the plot 4-3.

It is sufficient to establish two equations .

Let x₁ =0 the Ψ (0) = 1 = A + B and reading any other value from the graph ,for an arbitrary x₂ one has the second equation

Ψ (x₂) = A exp(ik₁ x₂ ) + B exp(-ik₁ x₂) . So A and B are found using the numerical values of Ψ(x). This numerical approach will prove to be useful when we examine the tunneling effect.

Results A(theoretical) = (1/2) ( 1 +i k₂ /k1) , B (theoretical) = (1/2) ( 1 -i k₂ /k1) .

Ath,Bth= 0.500E+00 0.500E+00 0.500E+00 -0.500E+00
A,B= 0.500E+00 0.509E+00 0.500E+00 -0.509E+00

real k1,k2
complex psi1, psi2 ,rooti ,A, B , Atheor, Btheor
c print*,'input : E, x1,x2 '
c read*,E , x1,x2
c print*,'psi1,psi2 '
c read*,psi1,psi2
data E,V,x1,x2/0.5,1.,0.,-1.05/
psi1=1.
psi2=1.38
k1=sqrt(2.*E)
k2=sqrt(2.*(V-E))
rooti=cmplx(0.,1.)
Atheor= (1./2.)*( 1. +rooti*k2/k1)
Btheor= (1./2.)*( 1. -rooti*k2/k1)
A=(psi1*exp(-rooti*k1*x2)-psi2*exp(-rooti*k1*x1))/
$(exp(rooti*k1*(x1-x2))-exp(rooti*k1*(x2-x1)))
B=(psi2*exp(rooti*k1*x1)-psi1*exp(rooti*k1*x2) )/
$(2.*rooti*sin(k1*(x1-x2)) )
print 100, atheor,btheor
print 110 ,a,b
100 format('Ath,Bth=',4(3x,e10.3))
110 format('A,B=',4(3x,e10.3))
stop
end

The case of E > V in fig 4-1 is considered next.

Let

k₁ = (2mE/h'² )^1/2 and k₂ = {2m(E-V)/h'² }^1/2 . (17)

In our units k₁ = (2E )^1/2 , k₂ = {2(E-V) }^1/2 .

We assume again that the original flux comes from the left toward the right hand side of the X axis. Ψ₁ (x) has the same form as before

Ψ₁ (x) = A exp( +ik₁x ) + B exp( -ik₁x ) . (18)

Ψ₂ (x) has the general form Ψ₂ (x) = C exp( +ik₂x ) + D exp( -ik₂x ) but we must set D=0 because it represents waves traveleing from

+ infinity to the left along the X axis. Accordingly

Ψ₂ (x) = C exp( +ik₂x ) . (19)

Just as before we have two continuity equations at the orgin (x=0) ,see eqs (3) -(4) , and three coefficients. Therefore, one of them is arbitrary.Let A=1 ,then the continuity equations give

1 + B = C , (20-a)

ik₁( 1-B) = ik₂ C . (20-b)

One gets

B = (k₁- k₂) /((k₁+ k₂) (21)

C= 2k₁/(k₁+ k₂ ) . (22)

We define the reflection coefficient ,R, as

R = B B^* /(A A^* ) = (k₁- k₂)² /((k₁+ k₂) ² . (23)

The transmission coefficient is T= CC^* /(A A^* ) = 1- R = { 2k₁/(k₁+ k₂ ) }² . (24)

Alternatively we can set C=1 in (19) and then the continuity conditions at x=0 give,

A + B = 1 , (25-a)

ik₁( A-B) = ik₂ , (25-b)

resulting in

A = (k₁ + k₂) /( 2k₁ ) , (26)

B = (k₁ -k₂)/( 2k₁ ) . (27)

Notice that A and B are real in the case E>V. Of course the same coefficients T and R result from (26 ) and (27). we have ,

T = CC^* /(A A^* ) = 1/A² = { 2k₁/(k₁+ k₂ ) }² , (28)

and R= B B^* /(A A^* ) = (k₁- k₂)² /(k₁+ k₂) ² . (29)

Since p = h'k , the coefficients T and R can be expressed in terms of linear momentum,

T = { 2p₁/(p₁+ p₂ ) }² , R = (p₁- p₂)² /(p₁+ p₂) ² .

Fig 4-4 Plot of Psi for E=2V=2 .

The current density is J = (ih'/(2m)){ Ψ ∂ Ψ^* /∂x - Ψ^* ∂ Ψ /∂x} = 1.41 , as seen from the numerical calculation given next.

x0,Real(j) Img(j)= -3. 1.41422522 -1.40756686E-008
x0,Real(j) Img(j)= 3. 1.41425657 1.78077215E-008

That is, the net flux on the left (x= -3) is equal to the flux on the right (x= +3.) ,namely 1.41 /time .

Fig 4-5 . Plot of Psi**2.

Since A and B are real when E > V we can relate the coefficients A and B to the maxima and minima of Ψ(x)² .

We have

Ψ(x)² = A² + B² + 2 AB cos(2k₁ x) . (30)

The maxima and minima are

Ψ(x)²_max. = (A+B)² , (31)

Ψ(x)²_min. = (A-B)² . (32)

Therefore

A =(1/2) { (Ψ(x)²_max. )^1/2 + (Ψ(x)²_min. )^1/2 } , (33)

and B =(1/2) { (Ψ(x)²_max. )^1/2 - (Ψ(x)²_min. )^1/2 } . (34)

Example: find A and B using Fig 4-5. We see that Ψ(x)²_max. =1 and Ψ(x)²_min.=0.5.

Eqs (33) and (34) give A =0.853 and B=0.146 . We corroborate these values using (26) and (27) .Given that E=2 and V=1 the values of k₁ and k₂ are , k₁ = 2 , k₂ = 2^1/2 and substitute them in (26)-( 27).

FORTRAN code

c plot of psi for potential barrier E > V

real k1,k2
complex psi,rooti ,aj
dimension Psi(0:5000)
data E,V/2.,1./
data xi,xf,nx/0.,6.,2000/
rooti=cmplx(0.,1.)
k1=sqrt(2.*E)
k2=sqrt(2.*(E-V))
dx=(xf-xi)/float(nx)
kp=int(float(nx)/40)
kount=kp
psi(0)=1.
psi(1)=psi(0)+dx*rooti*k2*psi(0)
print 100,xi,Real(psi(0)),aimag(psi(0))
do 10 i=2,nx
x=xi+dx*float(i)
if(x.le.0.)psi(i)=2.*psi(i-1)-psi(i-2) -dx**2*(2.*E)*psi(i-1)
if(x.gt.0.)psi(i)=2.*psi(i-1)-psi(i-2) -dx**2*(2.*(E-V))*psi(i-1)
if(i.eq.kount)then
print 100,x, real(psi(i)), aimag(psi(i))
kount=kount+kp
endif
10 continue
100 format('x,Re(psi), Im(psi)=',3(3x,e10.3))
print*,' '
x0=-3.
call density(x0,xi,dx,psi,raj,aimgaj)
print*,'x0,Real(j) Img(j)=',x0,raj,aimgaj
x0=+3.
call density(x0,xi,dx,psi,raj,aimgaj)
print*,'x0,Real(j) Img(j)=',x0,raj,aimgaj
stop
end

subroutine density(x0,xi,dx,psi,raj,aimgaj)
complex ajj , rooti,psi,dpsi
dimension Psi(0:5000)
rooti=cmplx(0.,1.)
i=int((x0-xi)/dx)
dpsi=(psi(i)-psi(i-1))/dx
raj=real( (rooti/2.)*( psi(i)*conjg(dpsi)-conjg(psi(i))*dpsi ) )
aimgaj=aimag((rooti/2.)*( psi(i)*conjg(dpsi)-conjg(psi(i))*dpsi) )
return
end

TUNNELING

Fig 4-6 . A flux of particles is incident from the left.

Particles with energy E < V have a probability of penetrating to the region x>a . This is called tunneling. One example of is the decay of alpha particles. Another is the conduction current between two metals separated by an insulator (a non conductor) ,see figure below.

Fig 4-7

The wave functions in the different regions are ( as before ,k= (2.*E)^1/2 , k₂= (2.*(V-E) )^1/2 ),

for x < 0 , Ψ ₁(x) = A exp(i kx) + B exp(-ikx) , (35)

for 0 < x < a , Ψ ₂(x) = Cexp(- k₂ x) + D exp(+k₂ x) , (36)

for x > a , Ψ ₃(x) = E exp(i k x) . (37)

we will set E=1. The constantsA, B , C and D are determined from the conditions of continuity of

Ψ and d Ψ /dx at the points x= 0 and a .

At x=a one has the linear system,

exp(- k₂ x) C + exp(+k₂ x) D = exp(i ka) (38)

k₂{- exp(- k₂ x) C + exp(+k₂ x) D } = ik exp(ika) (39)

Solving for c and D gives,

C == 1/2*(-I*k + k2)*e^(I*a*k + a*k2)/k2 , D == 1/2*(I*k + k2)*e^(I*a*k- a*k2)/k2 . (40)

SAGE code

C,D,k2,k,a=var('C,D,k2,k,a')
solve([e^(-k2*a)*C+e^(k2*a)*D==e^(I*k*a),-k2*e^(-k2*a)*C+k2*e^(k2*a)*D==I*k*e^(I*k*a)],C,D)

[[C == 1/2*(-I*k + k2)*e^(I*a*k + a*k2)/k2, D == 1/2*(I*k + k2)*e^(I*a*k
- a*k2)/k2]]

[[C == 1/2*(-I*k + k2)*e^(I*a*k + a*k2)/k2, D == 1/2*(I*k + k2)*e^(I*a*k - a*k2)/k2]]

At x=0 the continuity conditions yield

A + B = C+ D , (41)

ik(A-B) = k₂ ( -C+ D) . (42)

Solving for A in terms of C and D as given in (40) one obtains

A = { (k+ik₂)(-ik+k₂) exp(iak+ak₂) + (k-ik₂)(ik+k₂) exp(iak-ak₂) } /(4kk₂) . (43)

The transmission coefficient is a function of E and V and is given by

T = k E^* E /(k A^* A) = 1/ (A^* A) . (44)

The expression for T if manipulated algebraically gives ( see refs. 3 and 7)

T(E)=4.*(k*k2)**2/((k**2+k2**2)**2*sinh(a*k2)**2+4.*(k*k2)**2 ) . (45)

Fig 4-8. Transmission coefficient versus the ratio E/V , for a barrier of width a=1.

Fortran code

c tunneling coef
real k,k2
complex acoef,a1,a2,rooti
data a,Ei, Ef,ne,V/1. ,0.,1.,60,1./
D(e)=4.*(k*k2)**2/((k**2+k2**2)**2*sinh(a*k2)**2+4.*(k*k2)**2 )
rooti=cmplx(0.,1.)
de=(Ef-ei)/float(ne)
do 10 ie=0,ne
e=ei+de*float(ie)
k = (2.*E)**(1./2.)
k2 =(2.*(V-E))**(1./2.)
a1=(1./(4.*k2*k))*(k+rooti*k2)*(-rooti*k+k2)*exp(rooti*a*k+a*k2)
a2=(1./(4.*k2*k))*(k-rooti*k2)*(rooti*k+k2)*exp(rooti*a*k-a*k2)
acoef=a1+a2
print 100,E/V,D(E),1./abs(acoef)**2
10 continue
100 format('E/V,D,1./A**2=',3(3x,e10.3))
stop
end

Now we will cast the problem of the transmission coefficient as a numerical problem. Schrodinger's equation is solved

by the finite difference method where Ψ_n is obtained from the two previous values Ψ_n-1 and Ψ_n-2 ,

i.e. psi2=2.*psi1-psi0 -dx**2*{2.*(E-V)*psi1} ,

starting at x=a where Ψ(a) = E exp(ika) = exp(ika) and the derivative is dΨ(a) /dx = ik exp(ika) .

One integrates backwards towards x=0 . At that point A + B = Ψ(0) and ik(A-B) = dΨ(0) /dx , whence

A = (1/2) { Ψ(0) -(i/k) dΨ(0) /dx } (46)

Fig 4-9. Transmission coefficient by numerical integration of Schrodinger equation.

The pertinent FORTRAN code is given below.

c FORTRAN code for the numerical integration
c tunneling for E < V , backward integration
real k,k2
complex psi2,psi1,psi0,rooti,deriv ,acoef
data a,Ei,Ef,V/1.,0.,1.,1./
data ne,nx/40,1000/
rooti=cmplx(0.,1.)
xi=a
xf=0.
de=(ef-ei)/float(ne)
dx=(xf-xi)/float(nx)
kp=int(float(nx)/40)
kount=kp
do 20 ie=0,ne
e=ei +de*float(ie)
k=sqrt(2.*E)
k2=sqrt(2.*(V-E))
psi0=exp(rooti*k*a)
psi1=psi0+dx*rooti*k*psi0
do 10 i=2,nx
x=xi+dx*float(i)
psi2=2.*psi1-psi0 -dx**2*(2.*(E-V)*psi1)
deriv=(psi2-psi1)/dx
psi0=psi1
psi1=psi2
10 continue
acoef=(1./2.)*(psi2-(rooti/k)*deriv)
print 110,E/V,1./abs(acoef)**2
20 continue
100 format('x,psi=',2(3x,e10.3))
110 format('E/V,Tcoef=',2(3x,e10.3))
stop
end

Fig 4-10 . Ψ^* Ψ normalized such that Ψ (a) =exp(ika) , thus Ψ^* (a) Ψ (a) =1.

FORTRAN code

c FORTRAN code for the numerical integration
c tunneling for E < V , backward integration
real k,k2
complex psi2,psi1,psi0,rooti,deriv ,acoef
data a,Ei,Ef,V/1.,0.5,.5,1./
data ne,nx/40,1000/
rooti=cmplx(0.,1.)
xi=a
xf=-3.*a
de=(ef-ei)/float(ne)
dx=(xf-xi)/float(nx)
kp=int(float(nx)/40)
kount=kp
do 20 ie=1,1
e=ei +de*float(ie)
k=sqrt(2.*E)
k2=sqrt(2.*(V-E))
psi0=exp(rooti*k*a)
psi1=psi0+dx*rooti*k*psi0
print 100,a,abs(psi0)**2
do 10 i=2,nx
x=xi+dx*float(i)
if(x.lt.a.and.x.gt.0.)psi2=2.*psi1-psi0 -dx**2*(2.*(E-V)*psi1)
if(x.le.0.)psi2=2.*psi1-psi0 -dx**2*(2.*(E-0.*V)*psi1)
if(i.eq.kount)then
kount=kount+kp
print 100,x,abs(psi2)**2
print 100,a+Abs(a-x),1.
endif
deriv=(psi2-psi1)/dx
psi0=psi1
psi1=psi2
10 continue
acoef=(1./2.)*(psi2-(rooti/k)*deriv)
print 110,E/V,1./abs(acoef)**2
20 continue
100 format('x,psi**2=',2(3x,e10.3))
110 format('E/V,Tcoef=',2(3x,e10.3))
stop
end

Resonant tunneling

Fig 4-11 .The two barrier potential shows resonant tunnelig for E/V0 = 0.552

Fig 4-12 The transmission coefficient equals unity when E= 0.552V0.

Fortran code for resonant tunneling

c FORTRAN code for the numerical integration
c resonant tunneling for E < V , backward integration
real k,k2
complex psi2,psi1,psi0,rooti,deriv ,acoef
data a, Ei,V0/1.,0.01,2./
data ne,nx/60,1000/
Ef=V0
rooti=cmplx(0.,1.)
xi=3.*a
xf=0.
x1=a
x2=2.*a
x3=3.*a
de=(ef-ei)/float(ne)
dx=(xf-xi)/float(nx)
kp=int(float(nx)/40)
kount=kp
do 20 ie=0,ne
e=ei +de*float(ie)
k=sqrt(2.*E)
k2=sqrt(2.*(V0-E))
psi0=exp(rooti*k*xi)
psi1=psi0+dx*rooti*k*psi0
c print 100,a,abs(psi0)**2
do 10 i=2,nx
x=xi+dx*float(i)
psi2=2.*psi1-psi0 -dx**2*(2.*(E-V(x-dx,a,v0))*psi1)
c if(i.eq.kount)then
c kount=kount+kp
c print 100,x,abs(psi2)**2
c print 100,a+Abs(a-x),1.
c endif
deriv=(psi2-psi1)/dx
psi0=psi1
psi1=psi2
10 continue
acoef=(1./2.)*(psi2-(rooti/k)*deriv)
print 110,E/V0,log10(1./abs(acoef)**2)
20 continue
100 format('x,psi**2=',2(3x,e10.3))
110 format('E/V0,log10(T)=',2(3x,e10.3))
stop
end

function v(x,a,v0)
if(x.lt.0.)v=0.
if(x.ge.0..and.x.le.a)v=v0
if(x.gt.a.and.x.lt.2.*a)v=0.
if(x.ge.2.*a.and.x.le.3.*a)v=v0
if(x.gt.3.*a)v=0.
return
end

One dimensional scattering

An initial wave approaches from the right with amplitude exp(ikx). The energy is E= 1.2 and encounters a repulsive

barrier of height V=1 and width=a=1. Fig 4-13 shows that Ψ(x) is shifted forward relative to the original wave exp(ikx).

Fig 4-13 . With V=+1 , E=1.2 ,the scattered wave is shifted forward , in the region x > 0, relative to the plane wave exp(ikx).

FORTRAN code

c FORTRAN code for scattering probelm in one Dim by an attractive well
c V=+1 , E = 1.2
c solve numerically by backward integration
real k,k2
complex psi2,psi1,psi0,rooti,deriv ,acoef, psiplane
data a, V/1.,1./
data ne,nx/1,4000/
psiplane(x)=exp(rooti*k*x)
Ei=1.2
Ef=Ei
rooti=cmplx(0.,1.)
xi=-4.*a
xf= 5.*a
de=(ef-ei)/float(ne)
dx=(xf-xi)/float(nx)
kp=int(float(nx)/40)
kount=kp
do 20 ie=1,ne
e=ei +de*float(ie)
k=sqrt(2.*E)
k2=sqrt(2.*(E-V))
psi0=exp(rooti*k*xi)
psi1=psi0+dx*rooti*k*psi0
print 100,xi,real(psi0),real(psiplane(xi))
do 10 i=2,nx
x=xi+dx*float(i)
if(x.lt.0.)psi2=2.*psi1-psi0 -dx**2*(2.*(E-0.*V)*psi1)
if(x.ge.0..and.x.le.a)psi2=2.*psi1-psi0 -dx**2*(2.*(E-V)*psi1)
if(x.gt.a)psi2=2.*psi1-psi0 -dx**2*(2.*(E-0.*V)*psi1)
if(i.eq.kount)then
kount=kount+kp
print 100,x,real(psi2),real(psiplane(x))
c print 100,a+Abs(a-x),real(exp(rooti*k*(a+Abs(a-x))))
endif
deriv=(psi2-psi1)/dx
psi0=psi1
psi1=psi2
10 continue
c acoef=(1./2.)*(psi2-(rooti/k)*deriv)
c print 110,E/V,log10(1./abs(acoef)**2)
20 continue
100 format('x,Real(psi),Real(plane)=',3(3x,e10.3))
110 format('E/V,log10(T)=',2(3x,e10.3))
stop
end

Figure 4-14 shows that in the case of an attractive well (V = -1 , a=1), the wave is shifted backwards. It lags the original wave after scattering.