Bloch equation: adding gradients
I am not sure if this is related, but this cross-product is part of the differential equation, something like $\frac{\textrm{d}M}{\textrm{d}t} = \gamma M \times B$... For a true infinitesimal $\textrm{d}t$, this makes sense. For nummerical considerations, this needn't be true. Think of the following:


neglect relaxation for a moment
$\vec{M}_0$ has a certain magnitude (i.e. length)
when you apply an rf pulse, you expect a rotation of $\vec{M}_0$, and a vector's length is preserved under rotations
when you calculate the rotation from the cross product, with several time-steps $\Delta t$, you probably do something like $\vec{M}_{0,n+1} = \vec{M}_{0,n} + \Delta t\;\gamma\vec{M}_{0,n} \times \vec{B}$ (i.e. you compute the change from the previous magnetization and the rf pulse and add it to the previous magnetization vector)
this operation does not preserve the length of the vector: $\vec{M}_{0,n} \times \vec{B}$ is perpendicular to $\vec{M}_{0,n}$. When you add the vectors, the resulting vector is longer than $\vec{M}_{0,n}$ (Pythagorean theorem)


Hence, if you do not take appropriate countermeasures, your magnetization vector will grow. Maybe this is what you see in your simulation: At the end you have an insanely huge signal that does not belong there. You can easily check this via a simple test: Neglect relaxation and apply a huge flip angle, lets say $90000^\circ$ and monitor your magnetization vector length. Certainly it should remain constant. If it doesn't, there might be a bug in your code.

Also try to reformulate your calculations with rotations matrices and check if the results get better.