---

import numpy as np

def tepl_yav_sh(nx, nt, hx, ht, a, U_k):
    U = U_k.copy()
    for j in range(1, nt):
        for i in range(1, nx-1):
            U[i][j] = U_k[i][j-1] + ht*a/hx**2 * (
                U_k[i-1][j-1] - 2*U_k[i][j-1] + U_k[i+1][j-1]
            )
    return U

nx = 5
nt = 4
x1, xn = 0, 1
hx = (xn-x1)/(nx-1)
a = 1
ht = hx**2/2/a

X = [i*hx for i in range(nx)]
T = [j*ht for j in range(nt)]

U = [[1 for j in range(nt)] for i in range(nx)]
U = np.array(U, dtype=float)

U[0] = 24
U[nx-1] = 14
U[:, 0] = [24+66*x-76*x**2 for x in X]

U_k = U.copy()
c = 0
while c <= 4:
    c += 1
    U = tepl_yav_sh(nx, nt, hx, ht, a, U_k)
    U_k = U.copy()

print(U)