/*----------------------------------------*/
/*　4次ルンゲクッタ法の復習　　　　       */
/*----------------------------------------*/
#include<stdio.h>
#include <stdlib.h>
#include<math.h>

#define max_it     1000000     /*　最大時間ステップ数　*/
#define x_end      3e0         /*  計算終了のx */
void func(double x, double y, double *f);

main(void)
{
  int i,num;
  double x,y,f,k1,k2,k3,k4,f_fin,dx;
  FILE *output_file;
  
  output_file=fopen("output.dat","w"); /* 出力ファイルオープン*/ 

  /* ----------------  初期条件 ----------------- */

  x=0e0;
  dx=1e-2;
  y=0e0;

  /* -------------- 初期条件終わり ------------- */


  
  /* -------- 4次ルンゲクッタの時間積分 --------- */    

  for(num=0; num < max_it; num++){

    func(x,y,&f);
    k1=f*dx; /* 1段目 */

    func(x+dx*5e-1,y+k1*5e-1,&f);
    k2=f*dx; /* 2段目 */

    func(x+dx*5e-1,y+k2*5e-1,&f);
    k3=f*dx; /* 3段目 */

    func(x+dx,y+k3,&f);
    k4=f*dx; /* 4段目 */
      
    y = y + 1e0 / 6e0 * ( k1 + 2e0*k2 + 2e0*k3 + k4 );
    f_fin = 1e0 / 6e0 * ( k1 + 2e0*k2 + 2e0*k3 + k4 ) / dx;


    x=x+dx;
    
    /* --------            出力ルーチン              -----------*/

    if(num%3==0){   /* 3 step ごとに出力 */
      printf("number of iteration = %d\n",num);
  
      fprintf(output_file,"%f %f \n",(float)x,(float)y);
    }

    /* --------        出力ルーチン終わり              -----------*/

    /* --------        time step control ---------- */

    dx =1e-1;

    if(abs(f_fin) > 1e-5) {
      if( dx  > 1e-1 * fabs(y / f_fin)){
	dx =  1e-1 * fabs(y / f_fin) ;
      }
    }
    /* --------        time step control 終わり ---------- */

    if(x > x_end) break; /* 計算終了判定 */
  }

}

void func(double x, double y, double *f)
{

  *f=x*(y+1.);

}