u = zeros(Nx+3)
Ix = range(1, u.shape[0]-1)
# Boundary values: u[Ix[0]], u[Ix[-1]]
# Set initial conditions
for i in Ix:
u_1[i] = I(x[i-Ix[0]]) # Note i-Ix[0]
# Loop over all physical mesh points
for i in Ix:
u[i] = - u_2[i] + 2*u_1[i] + \
C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
# Update ghost values
i = Ix[0] # x=0 boundary
u[i-1] = u[i+1]
i = Ix[-1] # x=L boundary
u[i-1] = u[i+1]
Program: wave1D_dn0_ghost.py.