---

import numpy as np
import matplotlib.pyplot as plt

G = int(input("G = "))
S = int(input("S = "))
nx = int(input("nx = "))
nt = int(input("nt = "))

l = 1
alpha = 1

hx = l / (nx - 1)
X = [i*hx for i in range(nx)]

ht = hx**2 / (2*alpha)
T = [j*ht for j in range(nt)]

U = np.zeros((nx, nt))

fi_0 = G
fi_1 = S
psi = [G + (G + 3*S)*x - 2*(G + S)*x**2 for x in X]
F = 0

U[0] = fi_0
U[-1] = fi_1
U[:, 0] = psi

for j in range(1, nt):
    for i in range(1, nx-1):
        U[i][j] = (
            U[i][j-1]
            + ht*alpha/hx**2 * (
                U[i-1][j-1] - 2*U[i][j-1] + U[i+1][j-1]
            )
            + ht*F
        )

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

plt.figure(figsize=(10, 6))

for j in range(nt):
    if j == 0 or j == 1 or j % 10 == 0 or j == nt-1:
        plt.plot(X, U[:, j], label=f"t = {T[j]:.5f}", marker="o")

plt.xlabel("Длина сечения, м")
plt.ylabel("Температура, °C")
plt.title("Распределение температуры")
plt.legend(title="Время, с")
plt.grid(True)
plt.show()