1J = 1; % Value of the coupling J
2B = 0.1*J; % Value of the magnetic field
3Sx = [0 1;1 0]; % S_x operator
4Sz = [1 0; 0 -1]; % S_z operator
5I = eye(2); % Identity matrix
6Hspins = -J*kron(Sx,Sx)-B*(kron(Sx,I)+kron(I,Sx)); % Hamiltonian
7down = [0 1]'; % Quantum state down = [0 1]^T
8Psi_0 = kron(down,down); % Initial wavefunction
9T = 2*pi/(2*B); % Period of time
10Nt = 1000; % Number of steps to construct time vector
11ti = 0; % Initial time
12tf = 2*T; % Final time
13dt = (tf-ti)/(Nt-1); % Step time dt
14t = ti:dt:tf; % Time vector
15U = expm(-1i*Hspins*dt); % Time propagator operator U(dt)
16SSz = (kron(Sz,I)+kron(I,Sz))/2; % Operator S1^z + S2^z
17Mz = zeros(size(t)); % Average magnetization
18for n=1:length(t) % Iteration to find Psi(t) and Mz(t)
19 if n==1
20 Psi = Psi_0; % Initial wavefunction
21 else
22 Psi = U*Psi; % Wavefuntion at time t_n
23 end
24 Mz(n) = Psi'*SSz*Psi; % Average magnetization at time t_n
25end
26plot(t/T,real(Mz),'r-','LineWidth',3)
27xlabel('$t/T$','Interpreter','LaTex','Fontsize', 30)
28ylabel('$\langle M_z \rangle$','Interpreter','LaTex','Fontsize', 30)
29set(gca,'fontsize',21)