import numpy as np
from numpy.linalg import eig, solve

def matrix_power(A, n, x):
    eigenvalues, eigenvectors = eig(A)
    
    # Vector x will be a linear combination of
    # the eigenvectors: x = sum(a[i]@eigenvectors[i])
    a = solve(eigenvectors, x)
    
    size = len(eigenvalues)
    y = np.zeros(size)
    
    # Compute vector y = (A**n)@x
    for i in range(size):
        y = y + a[i]*(eigenvalues[i]**n)*eigenvectors[:, i]

    return y