ObviousNistrem4 (основной), неявный Runge Kutta 1 (вспомогательный) при переполнении стека

Вот что я сделал, он находит 3 первых значения с помощью одностадийного RK, но затем у меня есть несколько вопросов, как использовать Nistrema4th здесь.

В первой части я получаю eps а затем используйте одностадийный Runge Kutta и поместите эти 3 значения в массив как первое, второе и третье с оптимальным h (что является наименьшим шагом). Затем мне нужно начать с 4-го и использовать Nistrema, чтобы найти остальных.

double fx(double x, double y, double t)
{
return (2*x)-y+(t*t)-2*(sin(t)+1)+cos(t);
}
double fy(double x, double y, double t)
{
return x+(2*y)-sin(t)-(2*(t*t))+(2*t)-1;
}
double precise_x(double t)
{
return sin(t)+1;
}
double precise_y(double t)
{
return (t*t);
}

int main(int argc, char* argv[])
{
double h=0.025,h_opt, x[1000],y[1000],x_pre[10000],y_pre[10000],eps,current_x,current_y,x_prev,y_prev;
double a=0.,b=1.,t=a, x_prev1,y_prev1,current_x1,current_y1, h_temp[3],x1,y1;
int k=20,i;

cout << "Input value of the epsilon, please: ";
cin >> eps;
a = 0;  b = 1;  t = a; x[0] = 1; y[0] = 0;

current_x = x[0];
current_y = y[0];
int flag = 0;

current_x = x[0];
current_y = y[0];
while(!flag)
{

for(i=0;i<k;i++)
{
x_prev = current_x;
y_prev = current_y;

current_x = x_prev + (h*fx(x_prev,y_prev,t+h));
current_y = y_prev + (h*fy(x_prev,y_prev,t+h));

x[0]=x[1];
x[1]=x[2];
x[2]=x[3];
x[3]=current_x;

y[0]=y[1];
y[1]=y[2];
y[2]=y[3];
y[3]=current_y;
t+=h;
}
if((fabs(x_prev-current_x)+fabs(y_prev-current_y)<eps))
{
flag = 1;
}
else
{
h/=2;
t=0;
current_x=1.;
current_y=0.;
}
h_opt=h;
}
//===============???????????????????????????????????????????===
//use formulas of Nisterm (4th order ) to get the results
double n = 1./h_opt;
for(i=0;i<=n;i++)
{
x[i+4] = x[i+2] + (h_opt/3.0)*(8*fx(x[i+3],y[i+3],(i+3)*h)
-(5*fx(x[i+2],y[i+2],(i+2)*h))
+4*fx(x[i+1],y[i+1],(i+1)*h)
-fx(x[i],y[i],i*h));

y[i+4] = y[i+2] + (h_opt/3.0)*(8*fy(x[i+3],y[i+3],(i+3)*h)
-(5*fy(x[i+2],y[i+2],(i+2)*h))
+4*fy(x[i+1],y[i+1],(i+1)*h)
-fy(x[i],y[i],i*h));
}
t = a;
for(i=0;i<=n;i++)
{
cout<<"\n j "<<j;
x_pre[i] = precise_x(t);
y_pre[i] = precise_y(t);
t += h_opt;
}
cout << "x\ty\t\tprecise_x precise_y" << endl;
for(int i = 0; i <=k; i++)
{
cout  << x[i] << "\t" << y[i]<<"\t"<< x_pre[i] <<"\t"<<y_pre[i]<<endl;
}
system("PAUSE");
return 0;
}