def compile_higher_order_function(eqs, syms, params, order=2, funname='anonymous',
    return_code=False, compile=False):
    '''From a list of equations and variables, define a multivariate functions with higher order derivatives.'''

    from dolang.symbolic import stringify, stringify_symbol

    vars = [s[0] for s in syms]
    # TEMP: compatibility fix when eqs is an Odict:
    eqs = [eq for eq in eqs]

    if isinstance(eqs[0], str):
    # elif not isinstance(eqs[0], sympy.Basic):
    # assume we have ASTs
        eqs = list([ast.parse(eq).body[0] for eq in eqs])
        eqs_std = list( [stringify_symbol(eq, variables=vars) for eq in eqs] )
        eqs_sym = list( [ast_to_sympy(eq) for eq in eqs_std] )
    else:
        eqs_sym = eqs

    symsd = list( [stringify_symbol((a,b)) for a,b in syms] )
    paramsd = list( [stringify_symbol(a) for a in params] )
    D = higher_order_diff(eqs_sym, symsd, order=order)

    txt = """def {funname}(x, p, order=1):

    import numpy
    from numpy import log, exp, tan, sqrt
    from numpy import pi as pi_
    from numpy import inf as inf_
    from scipy.special import erfc

""".format(funname=funname)

    for i in range(len(syms)):
        txt += "    {} = x[{}]\n".format(symsd[i], i)

    txt += "\n"

    for i in range(len(params)):
        txt += "    {} = p[{}]\n".format(paramsd[i], i)

    txt += "\n    out = numpy.zeros({})".format(len(eqs))

    for i in range(len(eqs)):
        txt += "\n    out[{}] = {}".format(i, D[0][i])

    txt += """

    if order == 0:
        return out

"""
    if order >= 1:
        # Jacobian
        txt += "    out_1 = numpy.zeros(({},{}))\n".format(len(eqs), len(syms))

        for i in range(len(eqs)):
            for j in range(len(syms)):
                val = D[1][i,j]
                if val != 0:
                    txt += "    out_1[{},{}] = {}\n".format(i,j,D[1][i,j])

        txt += """

    if order == 1:
        return [out, out_1]

"""

    if order >= 2:
        # Hessian
        txt += "    out_2 = numpy.zeros(({},{},{}))\n".format(len(eqs), len(syms), len(syms))

        for n in range(len(eqs)):
            for i in range(len(syms)):
                for j in range(len(syms)):
                    val = D[2][n,i,j]
                    if val is not None:
                        if val != 0:
                            txt += "    out_2[{},{},{}] = {}\n".format(n,i,j,D[2][n,i,j])
                    else:
                        i1, j1 = sorted( (i,j) )
                        if D[2][n,i1,j1] != 0:
                            txt += "    out_2[{},{},{}] = out_2[{},{},{}]\n".format(n,i,j,n,i1,j1)

        txt += """

    if order == 2:
        return [out, out_1, out_2]

"""


    if order >= 3:
        # Hessian
        txt += "    out_3 = numpy.zeros(({},{},{},{}))\n".format(len(eqs), len(syms), len(syms), len(syms))

        for n in range(len(eqs)):
            for i in range(len(syms)):
                for j in range(len(syms)):
                    for k in range(len(syms)):
                        val = D[3][n,i,j,k]
                        if val is not None:
                            if val != 0:
                                txt += "    out_3[{},{},{},{}] = {}\n".format(n,i,j,k,D[3][n,i,j,k])
                        else:
                            i1, j1, k1 = sorted( (i,j,k) )
                            if D[3][n,i1,j1,k1] != 0:
                                txt += "    out_3[{},{},{},{}] = out_3[{},{},{},{}]\n".format(n,i,j,k,n,i1,j1,k1)

        txt += """

    if order == 3:
        return [out, out_1, out_2, out_3]
    """

    if return_code:
        return txt
    else:
        d = {}

        exec(txt, d)
        fun = d[funname]

        if compile:
            raise Exception("Not implemented.")

        return fun
From a list of equations and variables, define a multivariate functions with higher order derivatives.