#include <stdio.h>
#include <math.h>
#define IMAX 10000
double func(double x, double y);
void mk_plot_data(char *a, double x[], double y[], int n);
void mk_graph(char *f);

/*================================================================*/
/*       main function                                            */
/*================================================================*/
main(){
  double x[IMAX+10], y[IMAX+10];
  double final_x, h;
  double k1, k2, k3, k4;
  char temp;
  int ncal, i;


  /*--- set initial condition and cal range ---*/

  printf("\ninitial value x0 = ");
  scanf("%lf%c", &x[0], &temp);

  printf("\ninitial value y0 = ");
  scanf("%lf%c", &y[0], &temp);

  printf("\nfinal x = ");
  scanf("%lf%c", &final_x, &temp);

  do{
    printf("\nNumber of calculation steps( <= %d ) = ",IMAX);
    scanf("%d%c", &ncal, &temp);
  }while(ncal > IMAX);      


  /* --- size of calculation step --- */

  h=(final_x-x[0])/ncal;

     
  /* --- 4th Runge Kutta Calculation --- */

  for(i=0; i < ncal; i++){
    k1=h*func(x[i],y[i]);
    k2=h*func(x[i]+h/2.0, y[i]+k1/2.0);
    k3=h*func(x[i]+h/2.0, y[i]+k2/2.0);
    k4=h*func(x[i]+h, y[i]+k3);

    x[i+1]=x[i]+h;
    y[i+1]=y[i]+1.0/6.0*(k1+2.0*k2+2.0*k3+k4);
  }


  /* --- make a graph --- */

  mk_plot_data("out.txt", x, y, ncal);
  mk_graph("out.txt");

}


/*================================================================*/
/*       define function                                          */
/*================================================================*/
double func(double x, double y){
  double dydx;

  dydx=sin(x)*cos(x)-y*cos(x);

  return(dydx);
}


/*==========================================================*/
/*   make a data file                                       */
/*==========================================================*/
void mk_plot_data(char *a, double x[], double y[], int n){
  int i;
  FILE *out;

  out = fopen(a, "w");

  for(i=0; i<=n; i++){
    fprintf(out, "%e\t%e\n", x[i], y[i]);
  }

  fclose(out);
}

/*==========================================================*/
/*   make a graph                                           */
/*==========================================================*/
void mk_graph(char *f){
  FILE *gp;
     
  gp = popen("gnuplot -persist","w");
      
  fprintf(gp, "reset\n");
  fprintf(gp, "set terminal postscript eps color\n");
  fprintf(gp, "set output \"graph.eps\"\n");
  fprintf(gp, "set grid\n");


  /* -------  plat graph ---------*/

  fprintf(gp, "plot \"%s\" using 1:2 with line\n",f);
  fprintf(gp, "set terminal x11\n");
  fprintf(gp, "replot\n");
                    
  pclose(gp);
}