(* Mathematica program to evaluate the barrier potential in Quantum Mechanics. 
The barrier exists between x=-1 and x=1 with height V.
Are used units where \hbar=1, m=1, so that the Schroedinger equation reads

      -1/2 (d^2 u)/(d x^2) + ( V(x) -E ) u(x) =0

  date : 8/10/2009
  author: Jorge C. Romão
  email : jorge.romao@ist.utl.pt
*)

(* Definition of basic variables *)
k = Function[En, Sqrt[2 En]]
q = Function[{En, V}, If[En>V,Sqrt[2 (En-V)], I Sqrt[2 (V - En)]]]
w = Function[En, En]
pot=Function[{x,V},If[x<-1,0,If[x>1,0,V]]]

(* Definition of R, T, A and B *)
R=Function[{En,V}, I E^(-2 I k[En]) (q[En,V]^2-k[En]^2) Sin[2 q[En,V]]/
 (2 k[En] q[En,V] Cos[2 q[En,V]] - I (q[En,V]^2+k[En]^2) Sin[2 q[En,V] ])]
T=Function[{En,V},  E^(-2 I k[En]) 2 q[En,V] k[En]/
 (2 k[En] q[En,V] Cos[2 q[En,V]] - I (q[En,V]^2+k[En]^2) Sin[2 q[En,V] ])]
A=Function[{En,V},E^(I (q[En,V] -k[En])) ( T[En,V] E^(2 I (k[En] 
	+ q[En,V])) - E^(2 I k[En]) R[En,V] -1)/(-1 + E^(4 I q[En,V]))]

B=Function[{En,V},E^(I (q[En,V] -k[En])) ( R[En,V] E^(2 I (k[En] + q[En,V])) 
	- E^(2 I k[En]) T[En,V] + E^(2 I q[En,V]))/(-1 + E^(4 I q[En,V]))]

(* Definition of the wave function *)
psi=Function[{x,t,En,V},If[x<-1,(E^(I k[En] x) 
	+ R[En,V] E^(-I k[En] x)) E^(-I w[En] t), 
	If[x>1, T[En,V] E^(I k[En] x) E^(-I w[En] t),
	(A[En,V] E^(I q[En,V] x)+ B[En,V] E^(-I q[En,V] x) ) E^(-I w[En] t)]]]

(* Definition of the Real and Imaginary part of the wave function *)
Repsi=Function[{x,t,En,V},ComplexExpand[Re[psi[x,t,En,V]]]]
Impsi=Function[{x,t,En,V},ComplexExpand[Im[psi[x,t,En,V]]]]

(* Plots *)
(* Conventions Repsi[x, t, En, V], Impsi[x, t, En, V], pot[x,V] *)
Animate[Plot[{Repsi[x, t, 2, 1], Impsi[x, t, 2, 1],pot[x,1]}, 
            {x, -5, 5},PlotRange->{-2,2}],{t,0,100,0.1}]