#include <math.h>
#include <stdio.h>

double func1(double x, double y, double z);
double func2(double x, double y, double z);

void main()
{
        double *vx,*vy,sum,x,y,v,xx,yy,vv,tt,h,h2;
	double x0=0.,y0=0.7,v0=0.,rkm=20.;
	double f1,f2,f3,f4,g1,g2,g3,g4;
	float dt;
	int i,j,nstep;
	FILE *if1,*if2,*if3,*if4;
	char *fnam;
	scanf("%d",&nstep);
	scanf("%f",&dt);
	i=9000;
        vx= (double *) malloc(i*sizeof(double));
        vy= (double *) malloc(i*sizeof(double));
	fnam="Analytic";
	if((if1 = fopen(fnam,"w")) == NULL) printf("Error opening file_1 !!\n");
	fnam="Euler";
	if((if2 = fopen(fnam,"w")) == NULL) printf("Error opening file_2 !!\n");
	fnam="RK2";
	if((if3 = fopen(fnam,"w")) == NULL) printf("Error opening file_3 !!\n");
	fnam="RK4";
	if((if4 = fopen(fnam,"w")) == NULL) printf("Error opening file_4 !!\n");

	/* Analytic */
	for(i=0;i<nstep;i++){
          xx=i*dt;
          yy=y0*cos(sqrt(rkm)*xx);
          vv=-y0*sqrt(rkm)*sin(sqrt(rkm)*xx);
          vx[i]=xx;
          vy[i]=yy;
	  fprintf(if1,"%f %f\n",xx,yy);
	}

	/* Euler */
	h=dt;
	yy=y0;
	vv=v0;
	tt=x0;
	for(i=0;i<nstep;i++){
	  x=tt;
	  y=yy;
	  v=vv;
	  vv=vv+h*func1(x,y,v);
	  yy=yy+h*func2(x,y,v);
	  tt=tt+h;
	  vx[i]=tt-h;
	  vy[i]=yy;
	  fprintf(if2,"%f %f\n",tt-h,yy);
	}
	
	/* RK2 */
	h=dt;
	yy=y0;
	vv=v0;
	tt=x0;
	for(i=0;i<nstep;i++){
	  x=tt;
	  y=yy;
	  v=vv;
	  f1=func1(x,y,v);
	  g1=func2(x,y,v);
	  x=tt+h;
	  y=g1*h+yy;
	  v=f1*h+vv;
	  f2=func1(x,y,v);
	  g2=func2(x,y,v);
	  vv=vv+0.5*h*(f1+f2);
	  yy=yy+0.5*h*(g1+g2);
	  tt=tt+h;
	  vx[i]=tt-h;
	  vy[i]=yy;
	  fprintf(if3,"%f %f\n",tt-h,yy);
	}

	/* RK4 */
	h=dt;
	h2=0.5*h;
	yy=y0;
	vv=v0;
	tt=x0;
	for(i=0;i<nstep;i++){
	  x=tt;
	  y=yy;
	  v=vv;
	  f1=func1(x,y,v);
	  g1=func2(x,y,v);
	  x=tt+h2;
	  y=g1*h2+yy;
	  v=f1*h2+vv;
	  f2=func1(x,y,v);
	  g2=func2(x,y,v);
	  x=tt+h2;
	  y=g2*h2+yy;
	  v=f2*h2+vv;
	  f3=func1(x,y,v);
	  g3=func2(x,y,v);
	  x=tt+h;
	  y=g3*h+yy;
	  v=f3*h+vv;
	  f4=func1(x,y,v);
	  g4=func2(x,y,v);
	  vv=vv+0.166666666*h*(f1+2.*f2+2.*f3+f4);
	  yy=yy+0.166666666*h*(g1+2.*g2+2.*g3+g4);
	  tt=tt+h;
	  vx[i]=tt-h;
	  vy[i]=yy;
	  fprintf(if4,"%f %f\n",tt-h,yy);
	}
}
/*---------------------------------------------------*/	
double func1(double x, double y, double z)
{
	double rkm=20.,arg;
	arg=-rkm*y;
	return arg;
}
double func2(double x, double y, double z)
{
	double arg;
	arg=z;
	return arg;
}