Having issues towards the end where the functions to evaluate the symbolic matrices with given values at each time step of the ODE solver (using Euler method) being unable to properly substitute and evaluate the matrices with the given values.

# setup and solve equations
import numpy as np
from sympy import *
from sympy.physics.mechanics import *
from sympy.physics.mechanics import vlatex
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
init_vprinting()

# define our symbolic variables
m1, m2, m3, m4, m5, L1, L2, L3, L4, L5, g, t, m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, m41, m42, m43, m44 = symbols('m1 m2 m3 m4 m5 L1 L2 L3 L4 L5 g t m11 m12 m13 m14 m21 m22 m23 m24 m31 m32 m33 m34 m41 m42 m43 m44')

# optionally, set the masses and lengths
m1 = L1
m2 = L2
m3 = L3
m4 = L4

I1 = (m1*L1**2)/12
I2 = (m2*L2**2)/12
I3 = (m3*L3**2)/12
I4 = (m4*L4**2)/12

# dynamic symbols have implicit dependence on time
theta1, theta2, theta3, theta4, theta5 = dynamicsymbols('theta1 theta2 theta3 theta4 theta5')
theta1_dot = diff(theta1, t)
theta1_ddot = diff(theta1_dot, t)
theta2_dot = diff(theta2, t)
theta2_ddot = diff(theta2_dot, t)
theta3_dot = diff(theta3, t)
theta3_ddot = diff(theta3_dot, t)
theta4_dot = diff(theta4, t)
theta4_ddot = diff(theta4_dot, t)

theta = [theta1, theta2, theta3, theta4]
theta_dot = [theta1_dot, theta2_dot, theta3_dot, theta4_dot]
theta_ddot = [theta1_ddot, theta2_ddot, theta3_ddot, theta4_ddot]

# define coordinate transforms
Tsa = np.array([[cos(theta1), -sin(theta1), L1*sin(theta1)],
          [sin(theta1), cos(theta1), -L1*cos(theta1)],
          [0, 0, 1]])
Tab = np.array([[cos(theta2), -sin(theta2), L2*sin(theta2)],
          [sin(theta2), cos(theta2), -L2*cos(theta2)],
          [0, 0, 1]])
Tbc = np.array([[cos(theta3), -sin(theta3), L3/2*sin(theta2)],
          [sin(theta3), cos(theta3), -L3*cos(theta3)],
          [0, 0, 1]])
Tcd = np.array([[cos(theta4), -sin(theta4), L4*sin(theta4)],
          [sin(theta4), cos(theta4), -L4*cos(theta4)],
          [0, 0, 1]])


Tsb = Tsa @ Tab
Tsc = Tbc @ Tsb
Tsd = Tcd @ Tsc

# write p1, p2 in world frame (in center of rod)
p1 = Tsa @ np.array([0,L1/2,1])
p2 = Tsb @ np.array([0,L2/2,1])
p3 = Tsc @ np.array([0,L3/2,1])
p4 = Tsd @ np.array([0,L4/2,1])

# get the velocities
v1 = np.array([diff(p1[0],t), diff(p1[1],t)])
v2 = np.array([diff(p2[0],t), diff(p2[1],t)])
v3 = np.array([diff(p3[0],t), diff(p3[1],t)])
v4 = np.array([diff(p4[0],t), diff(p4[1],t)])

# total kinetic energy
# edit 2/15/24: fixed replaced 0.5*I2*theta2_dot**2 with 0.5*I2*(theta1_dot + theta2_dot)**2
T = m1*(v1.T @ v1)/2 + m2*(v2.T @ v2)/2 + m3*(v3.T @ v3)/2 + m4*(v4.T @ v4)/2 + (I1*theta1_dot**2)/2 + (I2*(theta1_dot + theta2_dot)**2)/2 + (I3*(theta1_dot + theta2_dot + theta3_dot)**2)/2 + (I4*(theta1_dot + theta2_dot + theta3_dot + theta4_dot)**2)/2

# total potential energy
V = m1*g*p1[1] + m2*g*p2[1] + m3*g*p3[1] + m4*g*p4[1]

# Lagrangian
L = T - V
# L = simplify(L) # failure to simplify often may result in no solution
L = expand(L)
L = trigsimp(L)
print(vlatex(L)) #print the equation 

# total kinetic energy
# edit 2/15/24: fixed replaced 0.5*I2*theta2_dot**2 with 0.5*I2*(theta1_dot + theta2_dot)**2
T = m1*(v1.T @ v1)/2 + m2*(v2.T @ v2)/2 + m3*(v3.T @ v3)/2 + m4*(v4.T @ v4)/2 + (I1*theta1_dot**2)/2 + (I2*(theta1_dot + theta2_dot)**2)/2 + (I3*(theta1_dot + theta2_dot + theta3_dot)**2)/2 + (I4*(theta1_dot + theta2_dot + theta3_dot + theta4_dot)**2)/2

# total potential energy
V = m1*g*p1[1] + m2*g*p2[1] + m3*g*p3[1] + m4*g*p4[1]

# Lagrangian
L = T - V
# L = simplify(L) # failure to simplify often may result in no solution
L = expand(L)
L = trigsimp(L)
print(vlatex(L)) #print the equation

row1_1 = collect(row1, theta1_dot)
row1_2 = collect(row1, theta2_dot)
row1_3 = collect(row1, theta3_dot)
row1_4 = collect(row1, theta4_dot)

row2_1 = collect(row2, theta1_dot)
row2_2 = collect(row2, theta2_dot)
row2_3 = collect(row2, theta3_dot)
row2_4 = collect(row2, theta4_dot)

row3_1 = collect(row3, theta1_dot)
row3_2 = collect(row3, theta2_dot)
row3_3 = collect(row3, theta3_dot)
row3_4 = collect(row3, theta4_dot)

row4_1 = collect(row4, theta1_dot)
row4_2 = collect(row4, theta2_dot)
row4_3 = collect(row4, theta3_dot)
row4_4 = collect(row4, theta4_dot)

print(vlatex(row1_1))
print(vlatex(row1_2))
print(vlatex(row1_3))
print(vlatex(row1_4))

print(vlatex(row2_1))
print(vlatex(row2_2))
print(vlatex(row2_3))
print(vlatex(row2_4))

print(vlatex(row3_1))
print(vlatex(row3_2))
print(vlatex(row3_3))
print(vlatex(row3_4))

print(vlatex(row4_1))
print(vlatex(row4_2))
print(vlatex(row4_3))
print(vlatex(row4_4))

M = Matrix([[row1_1, row1_2, row1_3, row1_4],[row2_1, row2_2, row2_3, row2_4],[row3_1, row3_2, row3_3, row3_4],[row4_1, row4_2, row4_3, row4_4]])
C = [[0, 0, 0, 0],
     [0, 0, 0, 0],
     [0, 0, 0, 0],
     [0, 0, 0, 0]]
for i in range(0, 4):
     for j in range(0, 4):
        for k in range (0,4):
            C[i][j] += (diff(M.row(i)[j], theta_dot[k]) + diff(M.row(i)[k], theta_dot[j]) + diff(M.row(k)[j], theta_dot[i]))*theta_dot[k]/2
for i in range(0,4):
     for j in range(0,4):
          print(vlatex(C[i][j])) 

B = [[1],
     [0],
     [0],
     [0]]
G = [0, 0, 0, 0]
for i in range(0, 4):
    G[i] = diff(V, theta[i])
print(vlatex(G))

B = Matrix(B)
G = Matrix(G)
C = Matrix(C)
M = Matrix(M)

# f_eqn = temp_inv*(B - C*theta_dot - G)
# f_eqn = f_eqn.subs([(theta1, 'theta1'), (theta2, 'theta2'), (theta3, 'theta3'), (theta3, 'theta3'), (theta4, 'theta4'), (theta1_dot, 'theta1_dot'), (theta2_dot, 'theta2_dot'), (theta3_dot, 'theta3_dot'), (theta4_dot, 'theta4_dot'), (cos, np.cos), (sin, np.sin)])

# print(shape(f_eqn))

def MassMf(y, d, j):
    y = Matrix(y)
    d = Matrix(d)
    # temp = M.evalf([(theta1, y[0]),(theta2, y[1]),(theta3, y[2]),(theta4, y[3]),(L1, d[0]),(L2, d[1]),(L3, d[2]),(L4, d[3]),(g,9.81),(t,j)])
    temp = M.evalf(subs = {theta1: y[0] , theta2: y[1], theta3: y[2], theta4: y[3], L1: d[0], L2: d[1], L3: d[2], L4: d[3], g: 9.81, t: j})
    return temp
def CoriolisMf(y1, y2, d, j):
    # print(y1)
    # print(y2)
    # print(d)
    # y1 = Matrix(y1)
    # y2 = Matrix(y2)
    # d = Matrix(d)
    # print(y1)
    # print(y1.row(0))
    # print("j: ")
    # print(j)
    # temp = C.evalf([(theta1, y1[0]),(theta2, y1[1]),(theta3, y1[2]),(theta4, y1[3]),(theta1_dot, y2[0]),(theta2_dot, y2[1]),(theta3_dot, y2[2]),(theta4_dot, y2[3]),(L1, d[0]),(L2, d[1]),(L3, d[2]),(L4, d[3]),(g,9.81),(t,j)])
    temp = C.evalf(subs = {theta1: y1[0], theta2: y1[1], theta3: y1[2], theta4: y1[3], theta1_dot: y2[0], theta2_dot: y2[1], theta3_dot: y2[2], theta4_dot: y2[3], L1: d[0], L2: d[1], L3: d[2], L4: d[3], g: 9.81, t: j})
    # print("temp: ")
    # print(temp)
    return temp
def GravityMf(y, d, j):
    # y = Matrix(y)
    # d = Matrix(d)
    # temp = G.evalf([(theta1, y[0]),(theta2, y[1]),(theta3, y[2]),(theta4, y[3]),(L1, d[0]),(L2, d[1]),(L3, d[2]),(L4, d[3]),(g,9.81),(t,j)])
    temp = G.evalf(subs = {theta1: y[0] , theta2: y[1], theta3: y[2], theta4: y[3], L1: d[0], L2: d[1], L3: d[2], L4: d[3], g: 9.81, t: j})
    return temp

# def ODE(t, y):
#     theta1 = y[0]
#     theta1_dot = y[1]
#     theta2 = y[2]
#     theta2_dot = y[3]

L = [1, 3, 4, 2]
tau = [0,0,0,0]
dt = 0.1
tf = 10
t = np.arange(0, tf, dt)
x0 = [0, np.pi/6, np.pi/4, np.pi/3, 0, 0, 0, 0]

# x_dot = [x2, (MassMf(theta,L)).inv()*tau - CoriolisMf(theta, theta_dot, L)*theta_dot - GravityMf(theta, L)]


x1 = [[0, np.pi/6, np.pi/4, np.pi/3]]
x2 = [[0,0,0,0]]

i = 0
j = 0

# print(x1.shape)
# print(x2.shape)

while j < tf:
    # print("x1: ")
    # print(x1)
    # print(x1[i])
    # print("x2: ")
    # print(x2)

    x1.append([x1[i][k] + dt * x2[i][k] for k in range(4)])
    # print(x1[i] + dt * x2[i])
    # print(x1)
    
    # x2.append(x2[i] + dt * (MassMf(theta,L)).inv()*tau - CoriolisMf(theta, theta_dot, L)*theta_dot - GravityMf(theta, L))
    # x2 = np.concatenate([x2, [np.transpose(x2[i]) + dt * ( - CoriolisMf(x1[i], x2[i], L) * np.transpose(x2[i]) - np.transpose(GravityMf(x1[i], L)))]], axis = 0)

    # print("x1: ")
    # print(x1)
    # print("x2: ")
    # print(x2)
    # print("j: ")
    # print(j)

    coriolis_matrix = CoriolisMf(x1[i], x2[i], L, j) #4x4
    gravity_vector = GravityMf(x1[i], L, j) #4x1
    # np_gravity_vector = np.array(gravity_vector)
    # coriolis_mult_result = np.array(coriolis_matrix) @ np.transpose(x2[i])
    # print(coriolis_mult_result)
    # coriolis_m_gravity = coriolis_mult_result - np_gravity_vector
    coriolis_m_gravity = coriolis_matrix * Matrix(x2[i]) - gravity_vector
    print(coriolis_m_gravity)
    x2.append([x2[i][k] + dt * coriolis_m_gravity[k] for k in range(4)])    

    i += 1
    j += dt
    print("END OF: " + str(i))
    

# f = [theta+dt*theta_dot, theta_dot + dt*x_dot[1]]

# sol = solve_ivp(x_dot, (0,tf), x0, t_eval=t)

Have tried using .subs function from SymPy with no success. Not sure on better method to evaluate all the values since each matrix cell is a very long equation.