/*

     解数值微分初值问题：

     龙格-库塔法求前k个初值 + 亚当姆斯法

 */

 #include<bits/stdc++.h>

 using namespace std;

 double f(double x,double y){

     //y(0) = 1

     return (y - *x/y);

 }

 void getRungeResult(double *Runge_k,double x0,double y0,double h,int N){

     //求解N个初值，保存在Runge_k[1 to N]中

     double K1,K2,K3,K4;

     double x1,y1;

     for(int i = ;i<=N;i++){

         x1 = x0+h;

         K1 = f(x0,y0);

         K2 = f(x0+h/,y0+h/*K1);

         K3 = f(x0+h/,y0+h/*K2);

         K4 = f(x1,y0+K4);

         y1 = y0 + h/*(K1+*K2+*K3+K4);

         Runge_k[i] = y1;

         x0 = x1;

         y0 = y1;

     }

     return;

 }

 //亚当姆斯多步法

 void Adams(double *Runge_k,double *predict,double x0,double y0,double h,int N){

     Runge_k[] = y0;

     //（0）龙格库塔法求前4个初值

     getRungeResult(Runge_k,x0,y0,h,);

     double y1,y2,y3,dy0,dy1,dy2,dy3;

     y1 = Runge_k[];

     y2 = Runge_k[];

     y3 = Runge_k[];

     dy0 = f(x0,y0);

     dy1 = f(x0+h,y1);

     dy2 = f(x0+*h,y2);

     dy3 = f(x0+*h,y3);

     double x3 = x0+*h;

     double x4,y4,yp,dyp,dy4;

     for(int i = ;i<=N;i++){

         x4 = x3+h;

         //(1)预测

         yp = y3 + h/*(*dy3-*dy2+*dy1-*dy0);

         predict[i] = yp;//保存预测值

         //预测要用dyp

         dyp = f(x4,yp);

         //(2)校正

         y4 = y3 + h/*(*dyp + *dy3 -*dy2+dy1);

         //存起来

         Runge_k[i] = y4;

         //求下一次需要用到导

         dy4 = f(x4,y4);

         //为下一次循环做准备

         x3 = x4;

         y3 = y4;

         dy0 = dy1;

         dy1 = dy2;

         dy2 = dy3;

         dy3 = dy4;

     }

     return;

 } 

 /*假设这里保证四阶精度*/

 int main(){

     /*说明：x0,y0是初值，h是小区间长度，N是要求的个数*/

     double x0,y0,h;

     int N;

     cout<<"输入初值x0,y0,小区间h,需要的初值个数N：";

     cin>>x0>>y0>>h>>N;

     //保存Runge求的4个初始值,龙格法求3个就可以;之后也用这个保存最终的Adams结果

     double Runge_k[];

     //保存预测值，方便以后比较

     double predict[];

     memset(predict,,sizeof(predict));

     memset(Runge_k,,sizeof(Runge_k));

     Adams(Runge_k,predict,x0,y0,h,N);

     cout<<endl;

     printf("预测值：");

     for(int i = ;i<=N;i++){

         if(i<){

             printf("%.6lf ",);

         }else{

             printf("%.6lf ",predict[i]);

         }

     }

     cout<<endl;

     printf("校正值：");

     for(int i = ;i<=N;i++){

         printf("%.6lf ",Runge_k[i]);

     }

 }