% Time base
dt=1e-4;steps=4e4;
t=[0:1:steps-1]*dt;
% Starting M
M1=[0,0,1]';
M0=[0,0,1]';
M=zeros(3,steps);M(:,1)=M1;
% Parameters of Experiment
B0=3;B1=1e-7;gamma=42.6e6;T1=2;T2=0.4;
% Set the pulse behavior
pulse=((t<0.5)|(t>2&t<2.5));
% This is the Larmor frequency
omega=2*pi*gamma*B0;
% Set the excitation frequency
omega0=omega+2*pi*2;
% dM=A M - M0 dt/T1
% omega0 is frequency of B1, omega is Larmor Frequency

% Time--step Matrix with field on
Aon=[-1/T2, (omega0-omega), 0;
    -(omega0-omega), -1/T2, 2*pi*gamma*B1;
   0, -2*pi*gamma*B1, -1/T1]*dt;

% With field off
Aoff=[-1/T2, (omega0-omega), 0;
    -(omega0-omega), -1/T2, 2*pi*gamma*0;
   0, -2*pi*gamma*0, -1/T1]*dt;

M0dtt1=M0*dt/T1;

% Generate the time history of M
for n=2:steps;
  if(pulse(n));A=Aon;else;A=Aoff;end;
  M(:,n)=A*M(:,n-1)+M0dtt1+M(:,n-1);
end;

% Plot the results
skip=10;
figure;plot(t(1:skip:end),M(:,1:skip:end)');grid on;
xlabel('t, time, sec.');ylabel('M, Net Magnetization');
legend('x','y','z');
