Let f(v(x)) be a composition of the functions f(v) and v(x), where both functions are multivariate. The first order chain rule applied on the composite function is simply the product of the Jacobians of the two functions, namely, fv*vx, where fv is the Jacobian of f and vx is the jacobian of v. This chain rule gives the first derivatives of f w.r.t x.
Higher order chain rules are much more complicated. In Levintal (2014) I derived compact notation for high order multivariate chain rules up to fifth order. These chain rules are implemented by the attached matlab package. The arguments of these chain rule functions are arrays of high order derivatives of f and v.

For example, chain3(fv,fvv,fvvv,vx,vxx,vxxx) calculates all third derivatives of f(v(x)) w.r.t x. f can be a vector-valued function. The arguments fv,fvv,fvvv are arrays that contain the first, second and third order derivatives of f w.r.t v, and vx,vxx,vxxx are similar arrays for v w.r.t x. For a fast way to obtain these derivatives see my function compderivs.m (which is more efficient than jacobian.m).

There are also optional arguments that can speed up the calculations. For more details, take a look at the readme file.

The package contains a test file that compares the chain rules to direct symbolic derivatives.

Reference:
Levintal, Oren, "Fifth Order Perturbation Solution to DSGE Models", 2014.