import numpy as np

# Coefficient matrix
A = np.array([[1, 2, 0, 0, 0],
              [2, -9, 3, 4, 0],
              [0, 3, -12, 4, 5],
              [0, 4, 4, -9, 5],
              [0, 0, 5, 5, -5]])

# RHS vector
b = np.array([1, 4, 9, 16, 25])

# Number of unknowns
n = 5

# Gaussian elimination method
for k in range(n-1):
    for i in range(k+1, n):
        factor = A[i, k]/A[k, k]
        b[i] = b[i] - factor*b[k]
        for j in range(k, n):
            A[i, j] = A[i, j] - factor*A[k, j]

# Back substitution
x = np.zeros(n)
x[n-1] = b[n-1]/A[n-1, n-1]
for i in range(n-2, -1, -1):
    sum = b[i]
    for j in range(i+1, n):
        sum = sum - A[i, j]*x[j]
    x[i] = sum/A[i, i]

print(x)