import scipy.optimize as resol
import numpy as np
import numpy.linalg as alg

A=1/9*np.array([[13,-16,-6],[-11,-1,3],[-13,-20,24]])
print( alg.eig(A))
B=1/25*np.array([[40,8,16],[40,-2,-4],[5,6,12]])

print( alg.eig(B))

def f(x):
    X=np.array(x)
    X=X.reshape((3,3))
    Y=A.dot(X)-X.dot(B)
    return Y.reshape(9)

X0=np.zeros(9)
X0[2]=1
sol=resol.root(f,X0)
Y=(sol.x).reshape((3,3))
print(Y)

X0=np.ones(9)
sol=resol.root(f,X0)
Y=(sol.x).reshape((3,3))
print(Y)