How to use the sympy.Function function in sympy

def g(v,w,i):
    from sympy import Lambda, Function ,symbols ,IndexedBase,Idx ,Max, Sum
    x = Function('x')
    i, n, j, dim, k =symbols('i, n, j, dim, k')
    v=IndexedBase('v')
    w=IndexedBase('w')
    net = Lambda((i, n, dim, k), Max(0.0, Sum(x(k)*w[n, k, i], (k, 0, dim-1))))
    dim = [10**4]*10
    index =[symbols('i%s'%m) for m in range(len(dim)+1)]
    y = [0]*len(dim)
    new = v[index[0]]
    for n in range(len(dim)):
        y[n] = net(index[n+1], n, dim[n], index[n])
        y[n] = y[n].subs(x(index[n]), new)
        new = y[n]
    return y[-1]

def function_lst(fstr, xtuple):
    sys.stderr.write(fstr + '\n')
    fct_lst = []
    for xstr in fstr.split():
        f = sympy.Function(xstr)(*xtuple)
        fct_lst.append(f)
    return(fct_lst)

def _divide_constriant(s, symbols, cons_index, cons_dict, cons_import):
    # Creates a CustomConstraint of the form `CustomConstraint(lambda a, x: FreeQ(a, x))`
    lambda_symbols = sorted(set(get_free_symbols(s, symbols, [])))
    r = generate_sympy_from_parsed(s)
    r = sympify(r, locals={"Or": Function("Or"), "And": Function("And"), "Not":Function("Not")})
    if r.has(Function('MatchQ')):
        match_res = set_matchq_in_constraint(s, cons_index)
        res = match_res[1]
        res += '\n        return {}'.format(rubi_printer(sympify(generate_sympy_from_parsed(match_res[0]), locals={"Or": Function("Or"), "And": Function("And"), "Not":Function("Not")}), sympy_integers = True))

    elif contains_diff_return_type(s):
        res = '        try:\n            return {}\n        except (TypeError, AttributeError):\n            return False'.format(rubi_printer(r, sympy_integers=True))
    else:
        res = '        return {}'.format(rubi_printer(r, sympy_integers=True))

    # First it checks if a constraint is already present in `cons_dict`, If yes, use it else create a new one.
    if not res in cons_dict.values():
        cons_index += 1
        cons = '\n    def cons_f{}({}):\n'.format(cons_index, ', '.join(lambda_symbols))
        if 'x' in lambda_symbols:
            cons += '        if isinstance(x, (int, Integer, float, Float)):\n            return False\n'
        cons += res

#
        # setup expression for the integration
        #

        step = Symbol('step')  # use symbolic stepsize for flexibility

        # let i represent an index of the grid array, and let A represent the
        # grid array.  Then we can approximate the integral by a sum over the
        # following expression (simplified rectangular rule, ignoring end point
        # corrections):

        expr = A[i]**2*psi_ho(A[i])*psi[i]*step

        if n == 0:
            print("Setting up binary integrators for the integral:")
            pprint(Integral(x**2*psi_ho(x)*Function('psi')(x), (x, 0, oo)))

        # But it needs to be an operation on indexed objects, so that the code
        # generators will recognize it correctly as an array.
        # expr = expr.subs(x, A[i])

        # Autowrap it.  For functions that take more than one argument, it is
        # a good idea to use the 'args' keyword so that you know the signature
        # of the wrapped function.  (The dimension m will be an optional
        # argument, but it must be present in the args list.)
        binary_integrator[n] = autowrap(expr, args=[A.label, psi.label, step, m])

        # Lets see how it converges with the grid dimension
        print("Checking convergence of integrator for n = %i" % n)
        for g in range(3, 8):
            grid, step = np.linspace(0, rmax, 2**g, retstep=True)
            print("grid dimension %5i, integral = %e" % (2**g,

# Need to replace all the strings correctly
            replicated_expression = replicated_expression.replace("%s{%d}" % (func.name, i + 1),
                                                                  "mod_solve[%d](x)" % i)

            # build function dict
            function_dict["%s_%d" % (func.name, i + 1)] = func

            # create sympy functions
            mod_solve.append(Function("%s_%d" % (func.name, i + 1)))

        # add the total flux at the end
        function_dict['total'] = composite_model

        replicated_expression += "- mod_solve[%d](x)" % num_models

        mod_solve.append(Function("total"))

        solutions = []
        # go through all models and solve for component fluxes algebraically
        for i, func in enumerate(composite_model.functions):
            solutions.append(solve(eval(replicated_expression), str(mod_solve[i]) + '(x)')[0])

        # use sympy to create new functions for the solved components
        component_flux = [lambdify(x, sol, function_dict) for sol in solutions]

        return component_flux

[-Derivative(u(t, x), t, t) + Derivative(u(t, x), x, x) == 0]

    References
    ==========

    .. [1] http://en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation

    """

    funcs = tuple(funcs) if iterable(funcs) else (funcs,)

    if not funcs:
        funcs = tuple(L.atoms(Function))
    else:
        for f in funcs:
            if not isinstance(f, Function):
                raise TypeError('Function expected, got: %s' % f)

    vars = tuple(vars) if iterable(vars) else (vars,)

    if not vars:
        vars = funcs[0].args
    else:
        vars = tuple(sympify(var) for var in vars)

    if not all(isinstance(v, Symbol) for v in vars):
        raise TypeError('Variables are not symbols, got %s' % vars)

    for f in funcs:
        if not vars == f.args:
            raise ValueError("Variables %s don't match args: %s" % (vars, f))

elif name == 'reduce':
                return Reduce(*args)

            elif name in _internal_zip_functions:
                return Zip(*args)

            elif name in _internal_product_functions:
                return Product(*args)

            else:
                msg = '{} not available'.format(name)
                raise NotImplementedError(msg)

        else:
            return Function(name)(*args)

    elif isinstance(expr, (int, float, Integer, Float, Symbol)):
        return expr

    else:
        raise TypeError('Not implemented for {}'.format(type(expr)))

except AttributeError:
            pass

        if len(exprs) != 4 or not all(len(exp.leaves) >= 1
                                      for exp in exprs[:3]):
            return

        if len(exprs[0].leaves) != len(exprs[2].leaves):
            return

        sym_args = [leaf.to_sympy() for leaf in exprs[0].leaves]
        if None in sym_args:
            return

        func = exprs[1].leaves[0]
        sym_func = sympy.Function(str(sympy_symbol_prefix + func.__str__()))(*sym_args)

        counts = [leaf.get_int_value() for leaf in exprs[2].leaves]
        if None in counts:
            return

        # sympy expects e.g. Derivative(f(x, y), x, 2, y, 5)
        sym_d_args = []
        for sym_arg, count in zip(sym_args, counts):
            sym_d_args.append(sym_arg)
            sym_d_args.append(count)

        try:
            return sympy.Derivative(sym_func, *sym_d_args)
        except ValueError:
            return

def replace_nosh_functions(expr):
    fks = []
    if isinstance(expr, float) or isinstance(expr, int):
        pass
    else:
        function_vars = []
        for f in expr.atoms(sympy.Function):
            if hasattr(f, 'nosh'):
                function_vars.append(f)

        for function_var in function_vars:
            # Replace all occurences of u(x) by u[k] (the value at the control
            # volume center)
            f = sympy.IndexedBase('%s' % function_var.func)
            k = sympy.Symbol('k')
            try:
                expr = expr.subs(function_var, f[k])
            except AttributeError:
                # 'int' object has no attribute 'subs'
                pass
            fks.append(f[k])

    return expr, fks

return_sensitivities (boolean): If True, return dictionary with
            sensitivity coefficients (as expressions).  Defaults to False.
        return_components (boolean): If True, return dictionary with
            component of evaluated uncertainty per argument/effect.
            Only considers uncorrelated part (which is the only thing
            implemented so far anyway).

    Returns:
        
        Expression indicating uncertainty
    """

    if collect_failures is None:
        collect_failures = set()
    import sympy
    u = sympy.Function("u")
    rv = sympy.sympify(0)
    sensitivities = {}
    components = {}
    for sym in recursive_args(expr):
        sym = aliases.get(sym, sym)
        try:
            sigma = sympy.diff(expr, sym)
            comp = sigma**2 * u(sym)**2
            rv += comp 
            sensitivities[sym] = sigma
            components[sym] = sympy.sqrt(comp)
        except ValueError as v:
            if on_failure == "raise" or not "derivative" in v.args[0]:
                raise
            else:
                warnings.warn("Unable to complete derivative on {!s}. "