如何用C++实现时滞方程的龙格－库塔法求解？

xsj3917

现有附件中的二维时滞方程，想用四阶龙格库塔法求解，如何用C++去实现？注意该方程的时滞项是变化的
附：一般微分方程的四阶龙格－库塔法C++代码：
for (i=0;i<n;i++)
{
x[i+1]=x+h;
k1=function(x,y);
k2=function(x+h/2,y+h*k1/2);
k3=function(x+h/2,y+h*k2/2);
k4=function(x+h,y+h*k3);
y[i+1]=y+h*(k1+2*k2+2*k3+k4)/6;
}
其中function为方程，例如

double function(double x,double y)
{
return y-2*x/y;
}

无水1324

建议用simulink进行仿真计算，比较简单

gghhjj

这个好像回答过，用隐式rk法，论坛好像能够找到相关程序