a     = 0.4; b     = 0.5; theta = 0.4 * pi/4;
dt    = 0.1; niter = 1000;

H = [ a^2*cos(theta)^2 + b^2*sin(theta)^2, (a^2-b^2)*cos(theta)*sin(theta);
	(a^2-b^2)*cos(theta)*sin(theta),      a^2*sin(theta)^2 + b^2*cos(theta)^2
	];
Udt = expm(-i * H * dt);
psit = [1.0; 0.0];
for iter = 1:niter
	t     = iter*dt;
	psit  = Udt * psit;
end
