Scattering solution to Schrödinger equation

Gaussian potential

Niels Walet, 1998-2009

Needs["PlotLegends`"]

Write the Schrödinger equation (mass is 1/2)

eqn=-(∂_(x,x))ϕ[x]+V[x] ϕ[x]-Energy ϕ[x]==0

Define the potential energy

V[x_]=-(E^(-(x^2)))

Scattering

Add these boundary conditions (at the left boundary to the equation). So we have an incoming plane wave. See that the wave gets a phase shift for positive x

scatter[Energy_,bnd_]:=Block[{eqa,solna},eqa={-(∂_(x,x))ϕ[x]+V[x] ϕ[x]-Energy ϕ[x]==0,ϕ[-bnd]==Sin[-(Sqrt[Energy]) bnd],ϕ'[-bnd]==(Sqrt[Energy]) Cos[-(Sqrt[Energy]) bnd]};solna=NDSolve[eqa,ϕ,{x,-bnd,bnd}];Plot[{V[x],ϕ[x]/.solna,Sin[(Sqrt[Energy]) x]},{x,-bnd,bnd},PlotLegend->{"V","ϕ","sine wave"}]]

Manipulate[scatter[Energy,3],{Energy,1,10}]

Bound States

We can also use a shooting method to find the bound state energy: if we start with a growing exponential on the left, we need a decaying exponential on the right

shoot[Energy_,bnd_]:=Block[{eqa,solna},eqa={-(∂_(x,x))ϕ[x]+V[x] ϕ[x]==Energy ϕ[x],ϕ[-bnd]==Evaluate[((Exp[(Sqrt[-Energy]) bnd])/25.)],ϕ'[-bnd]==Evaluate[(((Sqrt[-Energy]) Exp[(Sqrt[-Energy]) bnd])/25.)]};solna=NDSolve[eqa,ϕ,{x,-bnd,bnd}];
Plot[{V[x],ϕ[x]/.solna,ϕ'[x]/ϕ[x]/.solna},{x,-bnd,bnd},PlotLegend-> {"V","ϕ","ϕ'/ϕ"}]]

Find the energy that makes the logarithmic derivative flat

Manipulate[shoot[Energy,3],{Energy,-1,-0.001,Appearance-> Labeled}]