мааааааааааааленький код.py

import numpy
from numpy.linalg import solve

lambdas = [0.01, 0.05, 0.09] 
mu      = [1,    0.7,  1.5 ]
Ns      = [3,    3,    5   ]
N1s     = [2,    2,    3   ]

N = Ns[1]
N1 = N1s[1]
l1 = lambdas[2]

lams = [i * l1 for i in range(N, N-N1-1, -1)]
mus  = [(1+i)*mu[0] for i in range(N)]
print(mus)

lams = {k:i for k,i in enumerate(lams)}
mus  = {k:i for k,i in enumerate(mus)}

def coeff(i,j, lams, mus):
    if i == j:
        return -(lams.get(j, 0) + mus.get(j-1,0))
    if i == j+1:
        return mus.get(j,0)
    if i == j-1:
        return lams.get(i,0)
    else: return 0

mat = numpy.mat([
    [coeff(i,j, lams, mus) for i in range(N1+2)]
    for j in range(N1+2)
])

print(mat)
eigenvals, eigenvects = numpy.linalg.eig(mat)
idx = eigenvals.argsort()[::-1][0]
solnn = eigenvects[:,idx]
sol = solnn / numpy.sum(solnn)
print(sol)