I am new to coding in Python. My main task is to solve a PDE of the form
(https://i.sstatic.net/clwMaegY.png)
with periodic boundary conditions, say:
(https://i.sstatic.net/TpYx9QqJ.png)
so naturally i wanted to first solve the wave equation:
(https://i.sstatic.net/65JUhtOB.png)

The code i have come up with is this:

import matplotlib.pyplot as plt
import numpy as np
def initial_con(x):
    u = np.cos(2.* np.pi * x / L)
    return u
def initial_vel(x):
    return 0.

def solv(v, alpha):

u = np.zeros(((M + 1), (N + 1)))
f = np.zeros(N + 1)
g = np.zeros(N + 1)

for j in range (0, N+1): # --- Initial condition
    f[j] = initial_con(j * dx)
    g[j] = initial_vel(j * dx)
    u[0, j] = f[j]

for j in range (0, N+1): # --- First step ( u(-1)=u(1)-u'(0)*2*dt )
    if (j == 0):
        u[1, j] = (alpha ** 2) / 2 * (f[j+1] + f[N]) + (1. - alpha ** 2) * f[j] + dt * g[j]
    elif (j != 0 and j != N):
        u[1, j] = (alpha ** 2) / 2 * (f[j+1] + f[j-1]) + (1. - alpha ** 2) * (f[j]) + dt * g[j]
    else:
        u[1, j] = (alpha ** 2) / 2 * (f[0] + f[j-1]) + (1. - alpha ** 2) * f[j] + dt * g[j]

for i in range (2, M): # --- Time evolution
    for j in range (0, N+1):
        if j == 0:
            u[i+1, j] = (alpha ** 2) * (u[i, j+1] + u[i, N]) + 2. * (1. - alpha ** 2) * (u[i, j]) + u[i-1, j]
        elif (j != 0 and j != N):
            u[i+1, j] = (alpha ** 2) * (u[i, j+1] + u[i, j-1]) + 2. * (1. - alpha ** 2) * (u[i, j]) + u[i-1, j]
        else:
            u[i+1, j] = (alpha ** 2) * (u[i, 0] + u[i, j-1]) + 2. * (1. - alpha ** 2) * (u[i, j]) + u[i-1, j]

return u

"""Parameters."""
L = 1.                                   # --- Right boundary of the simulation domain
T = 2.                                   # --- Final time
M = 600                                  # --- Number of time steps
N = 200                                  # --- Number of space mesh points
v = 0.5                                  # --- Wave speed

dt = T / M
dx = L / N
alpha = v * dt / dx # --- CFL<1

u = solv(v, alpha)

x = np.linspace(0., L, N + 1)

plt.plot(x, u[0, :])
plt.ylabel("u(t0,x)")
plt.xlabel("x")
plt.show()

This code doesn’t seem to converge to the correct solution, which, given the initial conditions, should be something like:
(https://i.sstatic.net/46HVH9Lj.png)
I am not exactly sure where my mistakes are, however. Any suggestion would be helpful.