import sympy as sym

# Defining symbols:

x1, x2, u, a = sym.symbols('x1, x2, u, a')

# Defining the vectorial x:

x_vec = sym.Matrix([x1, x2])

# Defining the vectorial function:

f_vec = sym.Matrix([a*x1*x2, x1*u])

# Calculating the Jacobian:

J = f_vec.jacobian(x_vec)
print('J =', J)

# Calculating num value of J using lambdify():

fun_J = sym.lambdify([x1, x2, u, a], J, 'numpy')
val_J = fun_J(x1=1, x2=2, u=3, a=4)
print('val_J =', val_J)