// === Runge-Kutta法による微分方程式の初期値問題解析プログラム===
#include <stdio.h>
#include <math.h>
#define  EPS 1.0e-8
#define  N    50	// 分割数
#define  X0   0.0	// Xの初期値
#define  XN   5.0	// Xの終了値
#define  Y0   1.0	// Yの初期値

// 関数が主プログラムmainより後ろに書かれているときは、mainの前に
// このように定義を書いておく。そうしないと、関数が見つからないという
// エラーが表示される。
void S_RungK( double, double *, double, int, double(*)() );
double FUNC( double, double );

// 主プログラム
void main(){
   FILE *fout;
   int i;
   double x, h, y[101];

   fout=fopen("RK.csv","w");
   
   h=(XN-X0)/(double)N;
   x=X0;
   y[0]=Y0;
   printf(" X0=%f  XN=%f  Y0=%f  h=%f \n",x,XN,y[0],h);

   S_RungK(x,y,h,N,*FUNC);    
   for(i=0; i<=N; i++) {
      printf(" i=%f, y=%f \n", h*i,y[i]);
      fprintf(fout,"%f, %f\n", h*i,y[i]);
   }
   
   fclose(fout);
}

//** Function to calculate Runge-Kutta method **
//   入力; y[0]  =yの初期値
//         h     =xの刻み
//         N     =計算ステップ数（刻みｈでN回計算）
//         x     =xの初期値
//         FUNC  =dy/dxの関数
//   出力; y[N+1}=計算結果（N+1個の配列に入れられている）
void S_RungK( double x, double y[], double h, int n, double (*FUNC)()){
   int i; 
   double k1,k2,k3,k4;

   for(i=0; i<N; i++){
      k1=h*(*FUNC)(x, y[i]); 
      k2=h*(*FUNC)(x+0.5*h, y[i]+0.5*k1); 
      k3=h*(*FUNC)(x+0.5*h, y[i]+0.5*k2); 
      x=x+h;
      k4=h*(*FUNC)(x,y[i]+k3); 
      y[i+1]=y[i]+(k1+2.0*(k2+k3)+k4)/6.0;
   }
}

//    dy/dxの関数をここに書く
//    以下では、dy/dx=yの場合を与えている
double FUNC( double x, double y ){  
   return ( y );  
}