Использование dopri5 для построения системы ОДУ в матричной форме

Система уравнений, которые меня интересуют, выглядит следующим образом:

Я смог построить их, изменив пример, который кто-то опубликовал, выполнив следующие действия:

import scipy as sp
import pylab as plt
import numpy as np
import scipy.integrate as spi

#Constants
c13 = 4.2
c14 = 4.2
c21 = 4.3
c32 = 4.4
c34 = 4.4
c42 = 4.4
c43 = 4.4

e12 = 1.9
e23 = 2.5
e24 = 2.2
e31 = 2.0
e41 = 2.0

#Time
t_end = 700
t_start = 0
t_step = 1
t_interval = sp.arange(t_start, t_end, t_step)

#Initial Condition
r = [0.2,0.3,0.3,0.5]

def model(t,r):
    Eqs= np.zeros((4))
    Eqs[0] = (r[0]*(1-r[0]*r[0]-r[1]*r[1]-r[2]*r[2]-r[3]*r[3])-c21*((r[1]*r[1])*r[0])+e31*((r[2]*r[2])*r[0])+e41*((r[3]*r[3])*r[0]))
    Eqs[1] = (r[1]*(1-r[0]*r[0]-r[1]*r[1]-r[2]*r[2]-r[3]*r[3])+e12*((r[0]*r[0])*r[1])-c32*((r[2]*r[2])*r[1])-c42*((r[3]*r[3])*r[1]))
    Eqs[2] = (r[2]*(1-r[0]*r[0]-r[1]*r[1]-r[2]*r[2]-r[3]*r[3])-c13*((r[0]*r[0])*r[2])+e23*((r[1]*r[1])*r[2])-c43*((r[3]*r[3])*r[2]))
    Eqs[3] = (r[3]*(1-r[0]*r[0]-r[1]*r[1]-r[2]*r[2]-r[3]*r[3])-c14*((r[0]*r[0])*r[3])+e24*((r[1]*r[1])*r[3])-c34*((r[2]*r[2])*r[3]))
    return Eqs

ode =  spi.ode(model)

ode.set_integrator('dopri5')
ode.set_initial_value(r,t_start)
ts = []
ys = []

while ode.successful() and ode.t < t_end:
    ode.integrate(ode.t + t_step)
    ts.append(ode.t)
    ys.append(ode.y)

t = np.vstack(ts)
x1,x2,x3,x4 = np.vstack(ys).T

plt.subplot(1, 1, 1)
plt.plot(t, x1, 'r', label = 'x1')
plt.plot(t, x2, 'b', label = 'x2')
plt.plot(t, x3, 'g', label = 'x3')
plt.plot(t, x4, 'purple', label = 'x4')
plt.xlim([0,t_end])
plt.legend()
plt.ylim([-0.2,1.5])

plt.show()

Это, безусловно, дает мне сюжет, который я хочу. Однако я хочу завершить стохастический анализ с этим набором ОДУ, и по этой причине его гораздо проще смоделировать, если система ОДУ записана в матричной форме (таким образом, я могу легко изменить размерность ОДУ). шума и посмотреть, как это влияет на ОДУ). Я понимаю, как математически записать уравнение в матричной форме, но я не понимаю, как изменить свой код, чтобы в части «def model (t, r):» он читался как массив/матрица. Чтобы преобразовать уравнения в матричную форму, я могу определить:

b = np.array([1, 1, 1, 1]) 
A = np.array([[1, 1+c21, 1-e31, 1-e41], 
              [1-e12, 1, 1+c32, 1+c42],
              [c13+1, 1-e23, 1, 1+c43],
              [c14+1, 1-e24, 1+c34, 1]])

И тогда система уравнений будет (где x — вектор (x1,x2,x3,x4)):

х (т) = диаг (х) [b ^ {T}-Адиаг (х) х]

Итак, мой вопрос: как мне изменить место, где я определил свои ОДУ, чтобы я мог вводить их в виде матрицы вместо того, чтобы записывать каждое уравнение по отдельности? (это также упростит задачу, если я позже рассмотрю систему с более чем 4 измерениями)