Sympy trouble with simplifying partial derivatives

Question

I do the following to generate an expression in Sympy:

1. Create some matrix Q_{ij} which holds some functions \eta, \mu, \nu, of x and y.
2. Sum over indices and take some partial derivatives.
3. Substitute some simple expressions for \eta, \mu, and \nu (say \sin(x)*\cos(y)).
4. Try to simplify the expression so that it explicitly calculates the partial derivatives of the simple expressions.

When I do this, it gives me the following error: NotImplementedError: Improve MV Derivative support in collect

The specific code that I used is:

from sympy import *

# Set up system and generate functions
x, y, z = symbols('x y z')
i, j, k, m, p = symbols('i j k m p')
xi = Matrix([x, y, z])
lims = range(0, 3)

eta = Function('eta', real=True)(x, y)
mu = Function('mu', real=True)(x, y)
nu = Function('nu', real=True)(x, y)

Q = Matrix([[2/sqrt(3)*eta, nu, 0],
            [nu, -1/sqrt(3)*eta + mu, 0],
            [0, 0, -1/sqrt(3)*eta - mu]])

# Create complicated expression of partial derivatives
Phi_L1 = -sum(Eijk(3, p + 1, i + 1)
              *diff(diff(diff(Q[k, l], xi[j]), xi[j]), xi[p])
              *diff(Q[k, l], xi[i])
              for i in lims
              for j in lims
              for k in lims
              for l in lims
              for p in lims
              )

Phi_L1 = simplify(Phi_L1)

# Choose example functions and try to evaluate expression explicitly
eta1 = sin(x)*cos(y)
mu1 = sinh(x)*cosh(y)
nu1 = x**2*y**2
expr1 = Phi_L1.subs(eta, eta1).subs(mu, mu1).subs(nu, nu1)
simplify(expr1)

I couldn't find a simpler example that gives the same error. For example, the following works as intended:

f = Function('f', real=True)(x, y)
expr = diff(diff(f, x), y)
simplify(expr.subs(f, sinh(x)*cosh(y)))

Answer 1

'collect' and 'simplify' have known problems with higher order partial derivatives https://github.com/sympy/sympy/issues/9068 outlines the issue

The example shown is

import sympy as sp
x, y = sp.symbols("x y")
f, g = sp.symbols("f g", cls=sp.Function, args=(x,y))
f, g = f(), g()

expr1 = g.diff(x, 2) + f.diff(x, 2) + 5*f.diff(x, 2) + 2*f.diff(x) + 9*f.diff(x)
expr2 = g.diff(x, 2) + f.diff(x, 2) + 5*f.diff(x, 2) + 2*f.diff(x) + 9*f.diff(x) + g.diff(x).diff(y) + f.diff(x).diff(y) 

which works correctly and show the expected output for both expr1 and expo2

sp.collect(expr1, f) works wonderfully but sp.collect(expr2, f) fails with the known error as the implementation is not finished...