// wave.cpp: Wave motion on an elastic string
#include <iostream>
#include <math.h>
using namespace std;

const double v=25.;    // velocity of wave in m/s
const double L=1.;    // length in m
const double A=0.02;  // initial ampiltude in m

const int N=20; // number of spatial segments
const double dx=L/N;   // spatial discretization
const int kmax=80; // number of time steps to simulate
const double dt=0.001; // temporal discretization

double yk_1[N+1], yk[N+1], yk1[N+1]; // yi at time step k-1, k and k+1

int main(){

  for (int i=0; i<=N; ++i){   // all points
    double x = i*dx;
    yk_1[i]=A*sin(M_PI*x/L);  // initial displacment and boundary conditions
    yk[i]=yk_1[i];            // becasue initial velocity = 0
  }
  double a = pow( v*dt/dx, 2. ); 

  for (int k=1; k<=kmax; ++k){
    for (int i=0; i<=N; ++i) 
      cout << i*dx << " " << yk[i] << endl;     // output result
    cout << endl;

    // forward Euler
    for (int i=1; i<N; ++i)   // internal points
      yk1[i] = a * (yk[i+1] + yk[i-1]) + 2*(1-a)*yk[i] - yk_1[i]; 

    for (int i=1; i<N; ++i){
      yk_1[i] = yk[i];     // copy y(k) to y(k-1)
      yk[i] = yk1[i];      // copy y(k+1) to y(k), 
    }
  }
}