import sympy as sp

dX, dY, dZ, m, q0, q1, q2, q3 = sp.symbols('dX dY dZ m q0 q1 q2 q3')

f = sp.Matrix([dX + m * q0 + 1])

J = f.jacobian([dX, dY, q0])

print(J)