I tried to use a midpoint method and numerically solve the Schrödinger equation for the original Landau-Zener (LZ) problem: a $2\times 2$ Hamiltonian
$$\left(\begin{array}{c} \alpha t\\ \delta \end{array}\begin{array}{c} \delta \\-\alpha t \end{array}\right)$$
with initial condition $\psi=(1, 0)$ (the ground state) at some $t=−1000$, and say $\alpha=0.01$ and $\delta=0.04$. I took a huge amount of time slices $(10^8)$ which gives a time step of ~$10^{−5}$.
My goal is to reach the exact value from the LZ formula, but no matter how small a time step I take, I always have an error of 0.1%, after averaging over the oscillations that arise in the asymptotic behaviour.
Has anyone encountered this problem?
Here is my Matlab code:
alpha=0.01;
delta=0.04;
N=100000001;
ti=-1000;
tf=1000;
time=linspace(ti,tf,N);
sec=time(2)-time(1);
c1=[1,zeros(1,N-1)];
c2=[0,zeros(1,N-1)];
p=exp(-2*pi*delta^2/(alpha*2));
disp(sqrt(p))
for t=2:N
c1(t)=c1(t-1)-1i*sec*(alpha*(time(t-1)+sec/2)*(c1(t-1)-1i*sec/2*
(alpha*time(t-1)*c1(t-1)+delta*c2(t-1)))+delta*(c2(t-1)-1i*sec/2*(-alpha*time(t-
1)*c2(t-1)+delta*c1(t-1))));
c2(t)=c2(t-1)-1i*sec*(-alpha*(time(t-1)+sec/2)*(c2(t-1)-1i*sec/2*(-
alpha*time(t-1)*c2(t-1)+delta*c1(t-1)))+delta*(c1(t-1)-1i*sec/2*(alpha*time(t-
1)*c1(t-1)+delta*c2(t-1))));
end
disp(sum(abs(c1(end-(N-1)/50:end))/((N-1)/50+1)));
Also posted at Physics SE.