Не могу решить модель Джайлса-Атертона с Рунге-Куттом 4-го порядка; (

моя проблема заключается в следующем:
Мне нужно интегрировать Модель Джайлса-Атертона (описана в ответе)

Я использовал те же параметры, что и в ответе, но получил плохие результаты; (
Сначала в моей программе я увеличиваю H с 0 до 2500, затем уменьшаю с 2500 до -2500, а затем снова увеличиваю с -2500 до 2500.

Мой код:

#include <fstream>
#include <vector>
#include <cmath>
#include <algorithm>

using namespace std;

double dH = 0.05;
double OneSix = dH/6;

double fM(double M,double H,short int s,double* p){
double Ma = p[0] * (1/tanh((p[2]*M + H)/p[3]) - p[3]/(p[2]*M + H));
short int s2 = Ma - M < 0 ? -1 : 1;
double deltaM = 0.5*(1 + s2*s);
double dMadHe = p[0]*(p[3]/(p[2]*M + H)*(p[2]*M + H) -1/(p[3]*sinh((p[2]*M + H)/p[3])*sinh((p[2]*M + H)/p[3])));
return (deltaM*(Ma - M) + p[1]*p[4]*dMadHe*s)/(p[4]*s - p[2]*(deltaM*(Ma - M) + p[1]*p[4]*dMadHe*s));
}

void findloop(vector<double>& u,vector<double>& d, double* p){
double mu0 = 4 * 3.14 * 1e-7;
double Hend = 2500;
double H = 0;
double M = 0;
short int sign = 1;
for(int i = 1; i*dH < Hend; ++i){
H = i*dH;
double f1 = fM(M,H,sign,p);
double f2 = fM(M + f1*dH*0.5,H + dH*0.5,sign,p);
double f3 = fM(M + 0.5*dH*f2,H + dH*0.5,sign,p);
double f4 = fM(M + f3*dH,H + dH,sign,p);
M += OneSix * (f1 + 2*f2 + 2*f3 + f4);
}
sign = -1;
for(int i = 1; Hend - i*dH >= -Hend; ++i){
H = Hend - i*dH;
double f1 = fM(M,H,sign,p);
double f2 = fM(M + f1*dH*0.5,H + dH*0.5,sign,p);
double f3 = fM(M + 0.5*dH*f2,H + dH*0.5,sign,p);
double f4 = fM(M + f3*dH,H + dH,sign,p);
M += OneSix * (f1 + 2*f2 + 2*f3 + f4);
u.push_back(mu0*M);
}
reverse(u.begin(),u.end());
sign = 1;
for(int i = 1; i*dH - Hend <= Hend; ++i){
H = i*dH - Hend;
double f1 = fM(M,H,sign,p);
double f2 = fM(M + f1*dH*0.5,H + dH*0.5,sign,p);
double f3 = fM(M + 0.5*dH*f2,H + dH*0.5,sign,p);
double f4 = fM(M + f3*dH,H + dH,sign,p);
M += OneSix * (f1 + 2*f2 + 2*f3 + f4);
d.push_back(mu0*M);
}
return;
}

int main()
{
vector<double> up;
vector<double> down;
double params[5];
params[0] = 1.48 * 1e6;
params[1] = 0.0889;
params[2] = 0.000938;
params[3] = 470;
params[4] = 483;
findloop(up,down,params);
ofstream f1("up.txt");
ofstream f2("down.txt");
for(int i = 0; i < up.size(); ++i)
f1 << up[i] << endl;
for(int i = 0; i < down.size(); ++i)
f2 << down[i] << endl;
return 0;
}

В чем проблема?