線形系が縮退するアニメーションを死ぬほど汚いコードで書いた

やってること

大学の後輩が線形系の縮退のイメージが湧かない，と言っていたので2次元線形系の固有ベクトルが融合していくようなアニメーションを作ろうと思いました．解いている方程式は至極単純で，次の形をしています．

\frac{d}{dt}\Bigl(\begin{matrix}
x\\
y
\end{matrix} \Bigl)=
\Bigl(\begin{matrix}
-1 & -1\\
0 & -1-\epsilon
\end{matrix}\Bigl)\Bigl(\begin{matrix}x\\
y\end{matrix}\Bigl)

$2\times 2$行列の$(2,2)$成分だけに摂動$\epsilon$を入れておき縮退を予め解いてあります．そして，for文をぶん回して$\epsilon$を徐々に変化させ，そのステップ毎に複数の初期条件からscipyのodeintで積分していきます．
下のサンプルでは，安定なノードを縮退させています．もう少しfor文を回せば原点がサドルに変化するのも見れると思いますし，摂動の入れ方次第ではノードとスパイラルの転移も見れるはずです．

実際に書いたコード

反知性的なコードなので公開するか半日悩みましたが，誰かに有意義なコメントをしてもらうことに期待して公開します．for文の中でステップ毎に摂動パラメータ$e(=\epsilon)$をステップ数に応じて計算させ，それを元に微分方程式を解いています．

Linear_degenerate.ipynb

%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from numpy import sin, cos, pi
import matplotlib.animation as animation
# degeneration of linear system
def func(v, t, e):
    return [ -v[0] - v[1], -(1+e)*v[1] ] #  2-D linear system

# def func(v,t):
#    return [ -v[1] -0.5*v[0]*(v[0]**2+v[1]**2), v[0] -0.5*v[1]*(v[0]**2 + v[1]**2)]

# time length and unit step    
dt = 0.1
T = np.arange(0.0, 20, dt)
v0 = [[17,7], [17,17], [7,17], [-17,-7], [-17,-17], [-7,-17],[-7,17], [7,-17]]

# make animation
fig = plt.figure()

ims = []

for k in range(200):
    e = (150-k)/50
    #each initial state
    space = []
    for v, color in zip(v0, 'rbgrbgkk'):
        u = odeint(func, v, T, args=(e,))
        space += plt.plot(u[:,0], u[:,1], color)

    ims.append(space)
    
    
    '''u0 = odeint(func, v0[0], T, args=(e,))
    space0 = plt.plot(u0[:,0], u0[:,1], 'r')
    u1 = odeint(func, v0[1], T, args=(e,))
    space1 = plt.plot(u1[:,0], u1[:,1], 'b')
    u2 = odeint(func, v0[2], T, args=(e,))
    space2 = plt.plot(u2[:,0], u2[:,1], 'g')
    u3 = odeint(func, v0[3], T, args=(e,))
    space3 = plt.plot(u3[:,0], u3[:,1], 'r')
    u4 = odeint(func, v0[4], T, args=(e,))
    space4 = plt.plot(u4[:,0], u4[:,1], 'b')
    u5 = odeint(func, v0[5], T, args=(e,))
    space5 = plt.plot(u5[:,0], u5[:,1], 'g')
    u6 = odeint(func, v0[6], T, args=(e,))
    space6 = plt.plot(u6[:,0], u6[:,1], 'k')
    u7 = odeint(func, v0[7], T, args=(e,))
    space7 = plt.plot(u7[:,0], u7[:,1], 'k')

    ims.append(space0+space1+space2+space3+space4+space5+space6+space7)
    '''
plt.xlim(-17 ,17) 
plt.ylim(-17 ,17)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
ani = animation.ArtistAnimation(fig, ims, interval=27, blit=True, repeat=True)

plt.show()

実行例

$x$とは一次独立だった固有方向が段々$x$に潰れていくのが見えると思います．

問題点

先に記載したコードは問題だらけで，
1:for文中で同じような操作を何度も行っている
2:自動化したかったが初期値に応じた色の変更をどうしようか考えていなかった
3:せっかく縮退を見ているのに固有方向のプロットを忘れている
などが挙げられます（なら自分で直しておけ，という話ですが）．

もっと賢いコード が書けたらあとで何とかします． を教えてくれた優しい方のおかげで問題1,2が解決しました．改めて感謝します．解決していただいた問題の部分ですが，一応自分の書いた汚いコードもコメントアウトで残しておくことにしました．
初期条件の場所ももう少し工夫できたらな，と思いましたが系によって見やすい場所は違うので結局手で入れるしかないのかなぁと．