# -*- coding: utf-8 -*-
"""
Created on Thu Feb 23 13:37:48 2017

@author: lbnc
"""

import numpy as np
import matplotlib.pyplot as plt

N = 1000
nT = 400

c = 100.
dx = 0.2
dt = 0.9*dx/c
fac = c**2*dt**2/dx**2

xs = np.arange( N )*dx
ts = np.arange( N )*dt

x0 = 70
w0 = 10
w2 = 2*w0
y0 = np.exp( -( ( xs - x0 )/w2 )**2 )
ydot0 = (xs-x0)*c/(0.5*w2**2)*y0

ys = np.zeros( [N, nT ] )
ys[:,0] = y0 - dt*ydot0
ys[:,1] = y0

xi = np.arange( N-2 ) + 1
ti = np.arange( nT-2 ) + 1

for n in ti:
    
    for j in xi:
        
        ys[j,n+1] = 2*ys[j,n] - ys[j,n-1] + fac*( ys[j+1,n] - 2*ys[j,n] + ys[j-1,n] )

plt.plot(xs, ys[:,0], xs, ys[:,nT-1])
plt.ylabel("utslag")
plt.xlabel("posisjon")