Calculates the jacobian matrix for d/dx_i (f(x))

def jacobian(f, x, param):
    """
    Calculates the jacobian matrix for d/dx_i (f(x))
    """

    N = len(x)
    M = len(param)
    J = np.zeros((N, M))
    for i, x_i in enumerate(x):
        for j in range(M):
            parameters=[]
            for k, p_k in enumerate(param):
                    parameters.append(AdFloat(p_k, k == j))
            # parameters = [AdFloat(p_k, k == j) for k, p_k in enumerate(param)]
            val = f(AdFloat(x_i,0), parameters)
            J[i, j] = val.dx

    return J