import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve

def N(lmbd, mu):
    return lmbd/(mu - lmbd)

def U(lmbd, mu):
    return 1/(mu - lmbd)


L = 7
lambdas0 = [4, 10, 1, 3]
kappa = [1, 1, 2, 1, 4, 2, 3]
mu = [4, 9, 5, 2, 5, 1, 8]
Q = np.array([[0, 0.3, 0, 0.08, 0.02, 0.3, 0.2, 0],
              [0, 0, 0.5, 0, 0, 0.2, 0.3, 0],
              [0.04, 0.3, 0, 0.06, 0.2, 0.2, 0.2, 0],
              [0, 0.8, 0, 0, 0.01, 0, 0.14, 0.05],
              [0.4, 0.06, 0.02, 0.02, 0, 0.2, 0.3, 0],
              [0, 0, 0.1, 0.1, 0, 0, 0.5, 0.3],
              [0, 0, 0, 0.9, 0.1, 0, 0, 0],
              [0.06, 0, 0, 0.6, 0.14, 0, 0.1, 0.1]])

Omega = solve(Q - np.eye(L+1), 0)

plt.plot()