solving Bloch Eq. (ODE) with pumping

Hellow world!

I'm having some trouble with solving the Bloch Eq. when I add pumping to the Eqs.

The Bloch Eqs. without pumping term are defind as:

where M is the magnetization (in my case- of vapor) in the x/y/x axis, γ is the gyromagnetic ratio, and are relaxations rates.

I solve the Bloch Eq. using the following code:

gamma=350;
B=0.5;
T1=10;
T2=10;
tend=T1;
M0=[1;1;1]; %initial conditions
A=[-1/T2, gamma*B, 0; -gamma*B, -1/T2, 0; 0, 0, -1/T1]; %to represent the Bloch Eq. above using 
[t,M]=ode45(@(t, M) A*M [0 tend], M0);

The analitical solution in known and eqauls to

when FFT the results as follows

L=length(t);
Fs=1/mean(diff(t));
n=2^newtpow2(L);
Y=fft(M,n);
Y=abs(Y)
Y=Y(1:n/2m:);
f=(0:(n/2)-1)*Fs/n
f=f'
plot(f, Y(:,1));

I get a shift in the frequency for tend longer than ~5T_2

On the contrary, I get the correct frequency (with 2*pi factor in the x axis) for tend shorger tend (e.g. 0.1*T2), but the fft is distorted.

I appricaite help with this issue, but the following question concerns me even more

In the above Eqs. I assumed that the initial magnetizations in x/y/z are equal and the solution just describe its decay (and oscillations for Mx and My)

Now I want to add pumping which will add population to Mx and will oppose the relaxation in x.

I implement it by revising A:

A=[Rop-1/T2, gamma*B, 0; -gamma*B, -1/T2, 0; 0, 0, -1/T1]; 

But it seems to change My as well as Mx, and for high amplitude (e.g. Rop*10*T2) resutls seems to be very wierd.

Help please :)