Dt = (0.25*DeltaX)/c
# Dt = (16*pi)/3
rows = ceil((16*pi) / Dt)
cols = ceil((4*pi) / DeltaX + 3)

u = zeros((rows, cols))
v = zeros((rows, cols))

for j in range(cols):
    u[0][j] = cos(-DeltaX + j*DeltaX)
    v[0][j] = 0

for i in range(1, rows):
    for j in range(1, cols-1):

        # u(x,t): fill row (time) i of grid space for columns (space) j
        u[i][j] = u[i-1][j] + Dt*(v[i-1][j]) + (a/2)*(u[i-1][j+1] - 2*u[i-1][j] + u[i-1][j-1])

        # v(x,t): fill row (time) i of grid space for columns (space) j
        v[i][j] = v[i-1][j] + ((Dt*pow(c, 2))/pow(DeltaX, 2))*(u[i-1][j+1] - 2*u[i-1][j] + u[i-1][j-1])
        + (b/2)*(v[i-1][j+1] - 2*v[i-1][j] + v[i-1][j-1])

    u = apply_bcs(u, DeltaX, i)
    v = apply_bcs(v, DeltaX, i)