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

/****************** 図7.5の数値解を求めるプログラム *********************

	・波動方程式を解きます．

*************************************************************************/

#define Pi	4.*atan(1.)
#define eps	pow(10,-15)
#define C	1.
#define M	160
#define DT	0.00625
#define NT	10
#define NM	(int)(NT/DT)
#define XMIN	(double)0.
#define XMAX	(double)1.
#define DX	(double)((XMAX-XMIN)/M)

double u_t0(double x){
	return (cos(8.*Pi*x)+1.)/2.;
}

void solver(int n, double x[], double u[], double v[]){
	int j;
	if(n == 0){
		for(j=0; j<=M; j++){
			x[j] = XMIN+(double)j*DX;
			if(x[j] >= 3./8. && x[j] <= 5./8.)
				u[j] = u_t0(x[j]);
			else
				u[j] = 0.;
			v[j] = 0.;
		}
	}else{
		u[0] = 0.;
		u[M] = 0.;
		for(j=1; j<=M-1; j++){
			u[j] += DT*C*v[j];
		}
		for(j=1; j<=M-1; j++){
			v[j] += DT*C*(u[j-1]-2.*u[j]+u[j+1])/(DX*DX);
		}
	}
	return;
}

int main(void){

	int j, n=0;
	double t=0.;
	double x[M+1];
	double u[M+1];
	double v[M+1];

	while(n <= NM){
		solver(n, x, u, v);

		if(n%10 == 0)
			for(j=0; j<=M; j++)
				printf("%5.3lf, %12.7lf\n", x[j], u[j]);

		n++;
		t += DT;
	}

	return 0;

}