Unfortunately, chances are very low that you will actually find this. I will try to briefly explain why this is so.
Usually, what you want to achieve by fitting is to get an understanding of the parameters which make the M(H) curve look the way they do. This means that your model curve should have some physical relationship to the process of magnetization reversal going on in the specimen under investigation. [Otherwise you could simply use a set of suited [maybe orthogonal] functions and expand the M(H) curves in them, get a nice representation and learn nothing]. This means you need a physical model.
Now, in reality there are many different mechanisms by which magnetization reversal can occur. To mention a few: homogeneous rotation of magnetization (Stoner-Wohlfarth model), cration and (maybe hindered) motion of domain walls, curling, quantum tunneling of magnetization.... This alone should be reason enough to not think that there is anything like a "best equation".
Concerning ferro- and superparamagnetism, there is one more fundamental difference: ferromagnetism [in the field regime where there is hysteresis in M(H)] is inherently a non-equilibrium phenomenon [otherwise M(H) could not be 2- or multi valued, stability is just representing the fact that the barriers are so high, that thermodynamic equilibrium can only be reached on very long timescales, maybe exceeding the lifetime of our universe...]. In the discussion of superparamagnetism (no hysteresis), one assumes in modelling that equilibrium is given. In the presence of magnetic anisotropy, the correctness of the latter assumption may depend on the time scale involved with the measurement. What may look as an equilibrium M(H) curve of a superparamagnet with a slow probe (say, "ordinary" VSM or SQUID magnetometry) may actually reveal "blocking" for a fast probe of magnetic fluctuations (e.g. Mößbauer spectrometry).
Finally one last comment concerning superparamagnetism: Here, the "basic" curve is the Langevin function. The "danger" with this function is that many times, a superposition of only a few such curves produces a seemingly satisfactory representation of measured data. Nevertheless, in most cases one would still have no clue whether the parameters extracted from such a fit have anything to do with the physical reality in the specimen.
So - where to go from here? The answer is, that it may take MANY more experiments than just measuring some M(H) relationship in order to gain a sufficient understanding of magnetization reversal in a specimen. Once one has identified the mode of magnetization reversal, one may refine that understanding - oftentimes by even more experiments and additionally by modeling -- using a then appropriate model.
For ferromagnets there are interesting and intelligent measurement protocols which have been published over the years (one older author name I spontaneously remember is "Hesse", and there is quite a bit mof more recent work [watch out for terms like "Delta-m-method" and "FORC"]). Concerning superparamagnetism, you will find extremely useful references in the RG profiles of A. Tamion and F. Tournus.
I agree with the comment by Kai Fauth. There is no definite answer to this question. You should determine dominating magnetization reversal mode (if it exists) and try to use appropriate equation to describe M(H) hysteresis loop.
Dear Ahmed Hassan, I strongly disagree. The Langevin function is only best when it is adequate. The range of meaningful applicability is quite restricted. [The fact that many papers neglect this fact doesn't make this less true.]