f: single fluxion object or an array or list of fluxions, representing a scalar or vector function
        v: vector of variables in f with respect to which the Jacobian will be calculated
        v_mapping: dict mapping variables in f to scalar or vector of values

def jacobian(f, v, v_mapping):
    """
        f: single fluxion object or an array or list of fluxions, representing a scalar or vector function
        v: vector of variables in f with respect to which the Jacobian will be calculated
        v_mapping: dict mapping variables in f to scalar or vector of values

    """    
    # make sure f is in the form np.array([fl1, ...])
    if isinstance(f, Fluxion):
        f = [f]
    f = np.asarray(f)
    v = np.asarray(v)
    v_mapping = _check_input_vals(v_mapping)

    m = len(v)
    n = len(f)
    T = len(list(v_mapping.values())[0]) # number of values per variable
    
    J = np.zeros((m,n,T))
    for i, f_i in enumerate(f):
        seed = dict.fromkeys(v, 1)
            
        dfi_dvj = f_i.diff(v_mapping, seed)
        J[:,i,:] = dfi_dvj.T

    return J.squeeze()