I think that the physics of the problem are wrong. Unless I am mistaken, you are doing a phenomenological simulation of a partial hysteresis loop or IRM acquisition curve?
As it stands you are looking 5 (nrows) magnetic moments and assessing how their net magnetization (here the average moment) is progressively rotated to saturation with increasing field (B). The problem is that, physically, the net magnetization is not directly rotated as you do here. The individual magnetic moments are rotated and yield a resultant net magnetization.
What you need to do, for each grain with an initial angle w.r.t. to the field, rotate the moment according to the applied field. Then average all moments to get the net magnetization for that field step. Not that the average has to be a vector average, since the moments should themselves be vectors. To get a smooth curve, you will need to average over an ensemble of grains, 5 will not be enough. I would start with 100 to develop the code, but you will likely need >1000 to get a nice smooth curve.
Edit: I should also add, that you will need to make some assumption about your switching field (coercivity) distribution, otherwise if you are modeling hysteresis of IRM, then for identical grains it is just a step function.
I hope this gives you a good starting point to update your code.
