dts =  [ 0.032 0.016 0.008 0.004 0.002 0.001 0.0005 ]
norms = [];
clf
for dt = dts
% Create a grid to evaluate the basis function integrals in
n=10000;
x=linspace(-10,10, n);
h=x(2)-x(1);

% Length of supercell
L = 2;

% Initial k-point
k=0.1241;

% Exponent of wave function
a=3;

% Two reference basis functions, one even and one odd to mimic s and p.
wf = [ exp(-a*x.^2); x.*exp(-a*x.^2) ];

% Finite difference kinetic energy operator
e = ones(n,1);
T = -0.5*spdiags([e -2*e e], -1:1, n, n)/h^2;

% First index is the displacement (-L, 0, L), second is the 2x2=4 matrix
S = zeros(3,4); % Overlap, corresponds to St_MM
D = zeros(3,4); % Dipole, corresponds to Xt_MM
TT = zeros(3,4); % Kinetic energy
VV = zeros(3,4); % Potential energy
iter=0;
% Evaluate the matrix elements
for R=-1:1 % For displacelements (t in GPAW)
iter = iter + 1;
wf2 = [ exp(-a*(x-R*L).^2); (x-R*L).*exp(-a*(x-R*L).^2) ];
wfr2 = [ x.*exp(-a*(x-R*L).^2); x.*(x-R*L).*exp(-a*(x-R*L).^2) ];
wfT2 = [ exp(-a*(x-R*L).^2)*T  ; ((x-R*L).*exp(-a*(x-R*L).^2))*T ];
wfV2 = -[ cos(x*2*pi/L).*exp(-a*(x-R*L).^2); cos(x*2*pi/L).*(x-R*L).*exp(-a*(x-R*L).^2) ];
O = wf*wf2'*h;
PP = wf*wfr2'*h;
KK = wf*wfT2'*h;
VVtemp = wf*wfV2'*h;
D(iter, :) = PP(:);
S(iter, :) = O(:);
TT(iter, :) = KK(:);
VV(iter, :) = VVtemp(:);
end

% Calculate the initial matrix elements for ground state
X = 0;
Sk = zeros(2);
Tk = zeros(2);
Vk = zeros(2);
Dk = zeros(2);
for R=-1:1
    X = X + 1;
    f = exp(1j*L*R*k);
    Sk = Sk + reshape(S(X, :),2,2) * f;
    Tk = Tk + reshape(TT(X, :),2,2) * f;
    Vk = Vk + reshape(VV(X, :),2,2) * f;
    Dk = Dk + reshape(D(X, :),2,2) * f;
end

% Solve for the ground state
Hk = Tk + Vk
% Remove epsilon numerics to make perfectly Hermitian
Sk = (Sk' + Sk)/2;
Hk = (Hk' + Hk)/2;

[PSI, E] = eig(Hk,Sk)
C = PSI(:,1);
% C'*Sk*C = I
iter=0;

% Total simulation time
Ttot = 10;
% Total number of iterations
Ntot = round(Ttot / dt);

% Store wave functions along the propagation here
wf = zeros(2,Ntot);
% Store norm here
N = zeros(1,Ntot);

w=4.41; % Frequency of the laser pulse
A0=10; % Strength of electric field

k0=k; % Store the inital k-point
for iter = 1:Ntot
t = iter * dt;
% Laser induced vector potential at time t
A = A0*sin(t*w);
k = k0 + A;
dAdt = A0*w*cos(t*w); % Time derivative of vector potential

% Calculate new basis functions at updated k
Sk = zeros(2);
Tk = zeros(2);
Vk = zeros(2);
Dk = zeros(2);
X=0;
for R=-1:1
    X = X + 1;
    f = exp(1j*L*R*k);
    Sk = Sk + reshape(S(X, :),2,2) * f;
    Tk = Tk + reshape(TT(X, :),2,2) * f;
    Vk = Vk + reshape(VV(X, :),2,2) * f;
    Dk = Dk + reshape(D(X, :),2,2) * f;
end

Hk = Tk + Vk + dAdt*Dk;

% Propagator
G = expm(-1j * inv(Sk) * Hk * dt);

% Propagate wavefunctions
C = G*C;

% Store the norm and wavefunctions
N(:,iter) = C'*Sk*C;
wf(:, iter) = C;
end

% Plot the wavefunctions and norm
if 0
   subplot(2,1,1);
   hold off
   plot(real(wf)','r');
   hold on
   plot(imag(wf)','b');
   subplot(2,1,2);
   plot(real(N));
end

norms = [ norms std(real(N)) ];
hold off
loglog(dts(1:numel(norms)), norms,'o');
if numel(norms)>2
p = polyfit(log10(dts(2:numel(norms))), log10(norms(2:end)),1);
xlabel('log timestep');
ylabel('log std norm');
title(['Exponent ' num2str(p(1))]);
hold on
loglog(dts, 10.^polyval(p, log10(dts)),'-');
pause(0.01);
end
end

