N = 500; % number of simulated spins
M = 300; % length of a diffusion gradient
D = 2e-3; % "diffusion coefficient"

mx = ones(N,1); % Initial Magnetization points in x direction
my = zeros(N,1);
mz = zeros(N,1);
dp = (rand(N,2)-0.5)*1.8; % generate random spin distribution
dt = 4e-6; % time interval for bloch simulator
G = [4,0];  % set gradient amplitude

S = zeros(M,1);

%Simulate positive gradient
for n=1:M/2
    
    [mx,my,mz] = bloch(0,G,dt,800,800,0,dp,2,mx,my,mz); % simulate one time step
    dp = dp + randn(N,2)*D; % move spins randomly
    S(n) = mean(mx(:) + i*my(:)); % compute the FID
    figure(100), subplot(211),quiver(dp(:,1),dp(:,2),mx(:),my(:)); axis([-1,1,-1,1]);
    subplot(212), plot(abs(S));title('FID'), drawnow;
end

% Simulate negative gradient
for n=1:M/2
    
    [mx,my,mz] = bloch(0,-G,dt,800,800,0,dp,2,mx,my,mz); % simulate one time step
    dp = dp + randn(N,2)*D; % move spins randomly
    S(n+M/2) = mean(mx(:) + i*my(:)); % compute the FID
     figure(100),subplot(211), quiver(dp(:,1),dp(:,2),mx(:),my(:)); axis([-1,1,-1,1]);
   subplot(212), plot(abs(S));title('FID'); drawnow;  
end
