Passing matrices as input in scipy.integrate.solve_ivp python

I was solving a 2DOF spring-mass-damper system given below:

These are the 2 governing Equations

I have solved it in the following way:

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
m1 = 3
m2 = 5

k1 = 7
k2 = 9

c1 = 1
c2 = 2

f1 = 40
f2 = 4

start_time = 0
end_time = 60

initial_position_m_1 = 6
initial_velocity_m_1 = 0

initial_position_m_2 = 9
initial_velocity_m_2 = 4

delta_t = 0.1
def F(t, y):
    arr = np.array([
        y[1],
        (1/m1)*(f1*np.cos(3*t) - ((c1 + c1)*y[1] + (k1 + k2)*y[0]) + c2*y[3] + k2*y[2]),
        y[3],
        (1/m2)*(f2*np.sin(t**2) - c2*y[3] - k2*y[2] + c2*y[1] + k2*y[0])
    ])
    return arr

time_interval = np.array([start_time, end_time])
initial_conditions = np.array([initial_position_m_1, initial_velocity_m_1, initial_position_m_2, initial_velocity_m_2])

#######     solving the system of equations    ####
sol = solve_ivp(F, time_interval, initial_conditions, max_step = delta_t)

T = sol.t
Y = sol.y

Now, this is done by converting the 2 governing equations into 4 equations like this:

The problem with this is that I have to write each and every equation separately (as the function F)

Matlab has a way of solving it just with matrices using Ode45 function i.e. you don't have to write all the equations seperately in the function F in Matlab. You can enter the mass, stiffness and damping coefficients as matrices in it. Like this:

I am trying to solve a problem involving 30x30 matrices and if I do it in the above way, I will have to write 60 separate equations for the function F whereas in Matlab, I can pass the previously calculated 30x30 matrices directly into function. Is there any way of doing the same with solve_ivp in python or any such functions?

Thank you.

arr = np.array([
        y[1],
        (1/m1)*(f1*np.cos(3*t) - ((c1 + c1)*y[1] + (k1 + k2)*y[0]) + c2*y[3] + k2*y[2]),
        y[3],
        (1/m2)*(f2*np.sin(t**2) - c2*y[3] - k2*y[2] + c2*y[1] + k2*y[0])
    ])

can be written as (roughly):

f = np.array([0, f1*np.cos(3*t),0,f2*np.sin(t**2)])
M = np.array([
        [0, 1, 0, 0],
        [(k1+K2), (c1+c1), k2, c2],
        [0,0,0,1],
        [k2, c2, ....]])
arr = f[:,None] + M.dot(y)

That M array could be passed via args=(M,) (it's independent of t and y). Or just be a global to the function.