---

import numpy as np
import matplotlib.pyplot as plt

n = 5
x0, xn = 0, 1
dx = (xn - x0) / (n - 1)
X = [x0 + i * dx for i in range(n)]

P = [round(1 / (1 + 4*x*(1-x)), 4) for x in X[1:-1]]
P.insert(0, 0)
P.append(0)

F = [round(-2/25 * (2 - x*(1-x)/(1 + 4*x*(1-x))), 4)
     for x in X[1:-1]]
F.insert(0, 0)
F.append(0)

A = np.zeros((n, n))
A[0][0] = 1
A[n-1][n-1] = 1

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

Y = np.dot(np.linalg.inv(A), F)

print(np.round(Y, 4))

plt.plot(X, Y)
plt.show()