% propagation under a static Hamiltonian 

d=750;    % dipolar coupling
sw=1600;  % spectral width
sys=[1/2 1/2];  % two spin-1/2 spin system
tsteps=64; % points in FID

% Hamiltonian
H=d*(2*S(sys,1,'z')*S(sys,2,'z')-S(sys,1,'x')*S(sys,2,'x')-S(sys,1,'y')*S(sys,2,'y'));
sigma0=S(sys,1,'x');  % initial density matrix
detect=F(sys,'x');    % detection operator

dwell=1/sw;

sigma=sigma0; % initial density matrix
U=expm(-i*2*pi*dwell*H); % dwell time propagator
FID=zeros(1,tsteps);  % data array
for j=1:tsteps  % propagation loop
   FID(j)=trace(detect*sigma);  % detection
   sigma=U*sigma*U';  % propagation
end
FID1=FID;

% EXAMPLE 1

[V,evals]=eig(H); % eigenvectors and eigenvalues
phases=exp(-2*pi*i*diag(evals)*dwell); % convert eigenvalues to a vector in angular units
detectH=V'*detect*V; % detection operator in H eigenbasis
sigma0H=V'*sigma0*V; % initial density matrix in H eigenbasis
A= sigma0H.* (detectH.'); % amplitudes matrix
FID=zeros(1,tsteps);
n=length(evals); % dimensionality of Hilbert space
FID=FID+trace(A); % add constant terms
for j=1:n
   for k=j+1:n % loop over $k>j$ transitions
      v=2*A(j,k); % signal amplitude
      if abs(v)>1e-10 % ignore transitions of negligible amplitude
         p=phases(j)*conj(phases(k)); % phase factor
         for t=1:tsteps % loop over dwell times
            FID(t)=FID(t)+real(v); % accummulate FID
            v=v*p; % time evolution of signal
         end
      end
   end
end
FID2=FID;

tscale=[0:tsteps-1]*dwell*1000;
plot(tscale,real(FID1),'r-',tscale,real(FID2),'b*');
xlabel('time / ms');