Quantum Mechanics
Step potential with average energy
Simulation With GNU Octave

%*******************************************
% Wavepacket propagation using   step potential
% and exponential of H
% 
% U(x)  = the step potential with average energy : 
% 
%*******************************************

L = 100;
N = 400;
x = linspace(0,L,N)';
dx = x(2) - x(1);

% I do not have this in Octave (it should be in package ..)

function h = heaviside (x)
  h = zeros (size (x));
  h(x > 0) = 1.0;
  h(x == 0) = 0.5;
end



ko = 2;
a = 20;
dk = 2*pi/L;
km=N*dk;
k=linspace(0,+km,N)';

phi = exp(-a*(k-ko).^2).*exp(-i*6*k.^2);
psi = ifft(phi);
psi = psi/sqrt(psi'*psi*dx);

avgE = phi'*0.5*diag(k.^2,0)*phi*dk/(phi'*phi*dk);
U = avgE*heaviside(x-(L/2));
 
e = ones(N,1);
Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2;
H = -(1/2)*Lap + spdiags(U,0,N,N);
 
NT = 210; 
TF = 30;
T = linspace(0,TF,NT); 
dT = T(2)-T(1); 
hbar = 1;
 
E = expm(-i*full(H)*dT/hbar); 
 
 
 
%***************
%Simulate P(x,T) 
%***************


for t = 1:NT; 
 
psi = E*psi; 
P = conj(psi).*psi; 

% coloring black and red   ..

  if t < NT/2, plot(x,P,'k');
  else   plot(x,P,'r');
  end
 
filename=sprintf('output/%03d.png',t);
print(filename, "-dpng");

axis([0 L 0 0.15]);
xlabel('x (m)', 'FontSize', 16);
ylabel('probability density (1/m)','FontSize', 16);
pause(0.05);
end

%ffmpeg -i output/%03d.png heaviside.avi
 
% ------------------------------------ 
 
    The output:    
 heaviside.avi