Thanks for the interest, the full problem is the following. I have a differential equation defined in a rather strange way. Let n > 2, and denote the total differentiation with respect to ‘t’ by capital D. My differential equation is

D^(n+1) u(t) = f(t,u(t)) + D^n { g(t,u(t)) },

where u is a vector and t is a scalar, f and g are also vector valued, the same dimension as u.

My biggest problem is that I need to calculate the right-hand side, which involves calculating the n-th derivative of g(t,u(t)). For this I need to use the generalised chain rule. By rewriting the equation into first order form, the values of D^k u(t), 0 <= k <= n become state variables, so in essence I need to calculate various partial derivatives of g(t,u) and combine them with D^k u(t).

The first issue is that I really need the full tensor form of the derivative, because they are used as multi-linear forms, multiplied by various D^k u(t) values. As far as I can tell Taylor expansion assumes that all arguments of the multi linear forms are the same (hence you get a polynomial)

The second issue is that I could not see an example how to parameterise the function by those shifted values, such that the derivative functions do not need to be recompiled every time a different shift is used (which is every call).

I can now see that ForwardDiff will work (in theory) even if it is very slow to compile. I was also thinking of just using computer algebra to do the derivation, but the function ‘g’ is not a nice algebraic expression.

As you can see this is highly non-standard and at the moment I think Julia is my best bet.