<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:

def higher_order_diff(eqs, syms, order=2):
    '''Takes higher order derivatives of a list of equations w.r.t a list of paramters'''

    import numpy

    eqs = list([sympy.sympify(eq) for eq in eqs])
    syms = list([sympy.sympify(s) for s in syms])

    neq = len(eqs)
    p = len(syms)

    D = [numpy.array(eqs)]
    orders = []

    for i in range(1,order+1):

        par = D[i-1]
        mat = numpy.empty([neq] + [p]*i, dtype=object)   #.append( numpy.zeros(orders))
        for ind in non_decreasing_series(p,i):

            ind_parent = ind[:-1]
            k = ind[-1]

            for line in range(neq):

                ii = [line] + ind
                iid = [line] + ind_parent
                eeq = par[ tuple(iid) ]
                mat[tuple(ii)] = eeq.diff(syms[k])

        D.append(mat)

    return D