I have no memory of this but the fail whale seems to come from here:
which is trying to say that indexing two levels deep I[i][1] is too hard for the grad=Dual method to differentiate. (Probably it should be an error, not @debug.)
You can work around it by not making deeply nested structures, e.g. making a tuple of powers for x, and a separate one for y.
Maybe. I mean evaluating polynomials with all coefficients is always going to be more straightforward than evaluating them with some subset of the possible coefficients – that’s what you’re trying to do with this list of powers. If the only restriction is m+n<4 then at best this is an optimisation which may save you a factor of 2, compared to just keeping the zero coefficients. I suspect there are much larger factors around.
This is a harder problem than polynomials, and the design is likely to be quite different – involving some tuple of functions, perhaps.