---

import numpy as np

def m_j(n, h, R):
    A = [[0 for j in range(n)] for i in range(n)]

    for i in range(n):
        for j in range(n):
            if i == j:
                A[i][j] = 2 * R[i] / h**2
            elif i == j + 1 or i == j - 1:
                A[i][j] = -R[i] / h**2

    return A


# Задача Эйлера о потере устойчивости сжатого стержня
n = 5
x0 = 0
xn = 1
h = (xn - x0) / (n - 1)

X = [i * h for i in range(n)]
R = [19 * x * (1 - x) for x in X]  # (14 + 24) / 2 = 19
A = m_j(n, h, R)

E = 1e-10                        # повышенная точность
C = 0                            # счётчик итераций

v0 = np.ones(n)
v0[0], v0[n - 1] = 0, 0
lambda0 = np.sqrt(sum(v**2 for v in v0))

A = np.array(A, dtype=float)
A[0], A[n - 1] = 0, 0
A[:, 0], A[:, n - 1] = 0, 0
A[0][0], A[n - 1][n - 1] = 1, 1

invA = np.linalg.inv(A)
k = lambda0

while k > E:
    C += 1

    y1 = np.dot(1 / lambda0, v0)
    v1 = np.dot(invA, y1)
    lambda1 = np.sqrt(sum(v**2 for v in v1))

    k = abs(lambda1 - lambda0)
    v0 = v1
    lambda0 = lambda1

P = round(1 / lambda1, 9)
Y = v0 / lambda0

print("Pкр =", P)
print("Y =", np.round(Y, 9))
print("Количество итераций C =", C)