from pylab import *
npoints=100; #number of points in x is npoints+1
L=1.
c2=1.    #sound velocity squared
deltax=L/npoints   # step in x
deltat=0.01   # step in time
x=empty(npoints+1)   #coordinates
u=empty(npoints+1)   #currrent displacements
v=empty(npoints+1)    #current velocities
newu=empty(npoints+1)   #new displacements
newv=empty(npoints+1)   #new velocities
for i in range(npoints+1):
 	x[i]=i*deltax
	u[i]=sin(pi*x[i]/L)   #initial displacements
	v[i]=0.               #initialk velocities
for j in range(100000):   #make 100000 time steps
	newu[0]=0.        #left boundary condition for displacement
	newv[0]=0.        #left boundary condition for velocity
	newu[npoints]=0.  #right boundary condition for displacement
	newv[npoints]=0.  #right boundary condition for velocity
	for i in range(1,npoints):  #now new values for inside points
		newu[i]=u[i]+v[i]*deltat
		newv[i]=v[i]+c2*(u[i-1]-2.*u[i]+u[i+1])*deltat
	for i in range(1,npoints):  #make current values from new values
		u[i]=newu[i]
		v[i]=newv[i]
	if (j % 1000) ==0:
		print(u[50])   #prints th displacement of the midle point each 1000 time steps