выдает какую то ерунду

Если писать ерунду, то оно и выдаст ерунду

Уравнения там для верхней y_h и нижней y_l ветвей. И это обычные функциональные зависимости y=y(x), а не параметрические!

x(это напряженность поля H) = -H_max...H_max = deltaH*(i-[(N-1)/2]-1),

deltaH= 2*H_max/(N-1) i=1...N, N=2*k+1, k - натуральное число (а у Вас индекс складывается с величиной, имеющей размерность!!! это могло бы быть верным только в одном случае: если шаг по x равен единичному значению напряженности, что вряд ли адекватно условиям, плюс, у Вас половина точек является отрицательной областью одной ветви, а другая половина - положительной областью другой ветви)

y_h=A*atan((x/H_cs+1)*_tan_)

y_l=A*atan((x/H_cs-1)*_tan_),

где A=2*M_s/Pi, _tan_=tan(Pi*M_rs/M_r/2)

В этой модели беда с длиной "усов". Они больше не смыкаются в определенной точке, а асимптотически устремляются к едрене фене. Либо с этим придется мириться, либо одно из двух. :) В модели с параметрическим заданием функции другая проблема - искажение областей максимального изменения намагниченности при высоких степенях косинуса. Нет в этом мире совершенства, увы.

Запостили бы, что ле, скриншот с реальными данными. Желательно - только кривую и координатные оси. Если уж не можете сохранить в текстовом формате. Вряд ли Вы сами сможете выбрать наиболее адекватный вид функции, а я мог бы _быстро_ подобрать с использованием готовых "спецсредств".

ЗЫ. Комментарии мурра смело пропускайте мимо ушей. Он несет отчаянную ахинею насчет того, что и как я понимаю. Да и вообще, не особо въезжает в тему. У него до сих пор все еще 100x100...