---

import numpy as np

def Пуассон(nx, nt, hx, ht, F, U_k):
    U = U_k.copy()
    for i in range(1, nx - 1):
        for j in range(1, nt - 1):
            U[i, j] = ((U_k[i-1, j] + U_k[i+1, j]) / hx**2
                       + (U_k[i, j-1] + U_k[i, j+1]) / ht**2
                       - F) / (2 / hx**2 + 2 / ht**2)
    return U

nx, nt = 5, 3
hx, ht = 1/4, 1/2
x1, xn = 0, 1
t1, tn = 0, 1

U = [[1 for j in range(nt)] for i in range(nx)]

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

fi = [96*t*(1-t) for t in T]
psi = [56*x*(1-x) for x in X]
F = -304

U = np.array(U, dtype=float)

U[0], U[nx - 1] = fi, fi
U[:, 0], U[:, nt - 1] = psi, psi

e = 0.1
z = np.max(np.abs(U))
U_k = U.copy()
k = 0

while z >= e:
    U = Пуассон(nx, nt, hx, ht, F, U_k)
    z = np.max(np.abs(U_k - U))
    U_k = U.copy()
    k += 1

k, z, np.round(U[::-1], 3)