How to convert symbolic expressions to an ODE and solve it in Python?

Question

Please quickly go through the question I asked in MATLAB Ask:
https://www.mathworks.com/matlabcentral/answers/827750-convert-the-symbolic-expression-into-the-expression-which-can-be-put-on-the-rhs-of-an-ode?s_tid=srchtitle
Now I want to reproduce the same work in Python and meet some problems. Given 2 state equations:

f1 = x2 
f2 = -c*sinx1-(a + b*cosx1)*x2

The sensitivity equation is given by

dSdt = PfPx*S+PfPw=RHS # dSdt: time derivative of S (sensitivity matrix),in this case,2x3 matrix 
PfPx # Jacobian matrix of F = (f1,f2) wrt x = (x1,x2),2x2 matrix 
S # Sensitivity matrix S = [[x3, x5, x7],[x4, x6, x8]],2x3 matrix 
PfPw # Jacobian matrix of F wrt parameters w = (a,b,c),2x3 matrix

By symbolic matrix calculation, RHS = PfPx*S + PfPw
Sensitivity equation is given as dSdt(i,j) = RHS(i,j) by vector form for dx3dt~dx8dt appending with state equation for dx1dt and dx2dt (see the picture at the end in that Link), when a=1,b=0,c=1, it generates 8 ODEs:

dx1dt= x2;
dx2dt = -sinx1-x2; # dx1dt~dx2dt are state eqns 
dx3dt = x4;
dx4dt = -x3cosx1-x2-x4;
dx5dt = x6;
dx6dt = -x5cosx1-x6-x2cosx1;
dx7dt = x8; 
dx8dt = -cos(x1)x7-x8-sinx1; # dx3dt – dx8dt are sensitivity eqns 

Try to code the symbolic calculation in Python,

import sympy as sym
import numpy as np
from scipy.integrate import odeint
    
def model(x, t): 
    a, b, c = sym.symbols('a,b,c')
    f1 = x[1]
    f2 = -c * sym.sin(x[0]) - (a + b * sym.cos(x[0])) * x[1]
    funcs = sym.Matrix([f1, f2])
    pfpx = funcs.jacobian([x[0], x[1]])  # pfpx = Jacobian('x1,x2,x3',[f1, f2, f3])
    pfpw = funcs.jacobian([a, b, c])
    s = sym.Matrix([[x[2], x[4], x[6]], [x[3], x[5], x[7]]])
    pfpx = pfpx.subs([(a, 1), (b, 0), (c, 1)])
    pfpw = pfpw.subs([(a, 1), (b, 0), (c, 1)])
    res = pfpx * s + pfpw
    dydt = [f1, f2, res[0], res[3], res[1], res[4], res[2], res[5]]
    return dydt

    tspan = np.linspace(0, 10)
    z0 = np.zeros(8)
    z0[0:2] = [1, 1]
    zrec = odeint(model, z0, tspan)

it is saying error on the last line zrec = odeint(model, z0, tspan) and the error msg is:

ValueError: 
Can't calculate derivative wrt 1.00000000000000. 

Any suggestions are appreciated!