Finding expression of $n$-th derivative, when $n$ is large For completeness, assume $C$ is an Archimedean copula with some generator function $\varphi$, which is usually assumed to have nice properties. It is known that $$ C(u_1, u_2, \ldots, u_n)=\varphi^{-1}(\varphi(u_1) + \ldots +\varphi(u_n))$$If my math is correct, for an $n-$variable copula, we can estimate its density function by 
$$ \frac{\partial^n C}{\partial u_1 \partial u_2 \cdots \partial u_n } = \{\varphi^{-1}\}^{(n)}(\varphi(u_1) + \ldots +\varphi(u_n))\prod_{j=1,\ldots,n} \varphi'(u_j) $$
However, with large $n$, it is not a surprise that it's quite computationally/memory intensive finding $\varphi^{-1}$ derivatives of the $n-$th order. 
My question: what are the most well known/fastest strategies for quick derivations of $\{\varphi^{-1}\}^{(n)}$? I'm not entirely sure what exactly I'm looking for here -- either a faster/cleaner way of finding the expression for density that wouldn't involve derivatives (if such exists?), or some sort of an approximation, that would make finding the expressions a more feasible task?
As a side note, I'm trying to get the expression of the density function and plug-in various values to it. Currently I've tried using deriv function on R, which is able to handle the derivations up to say $n \sim 15$. 
A possible workaraound would be to simulate a sample from the wanted copula distribution and non-parametrically estimate its density via some kernel smoother, but as far as I know, kernel smoothers also are quite complex and slow in large dimensions. 
 A: I am not an expert in this area, so please take this answer with a grain-of-salt.  Since you have not specified a functional form for your objects I am going to give an answer that applies at a broad level of generality.  (No doubt there are better methods if you are willing to specify more about your particular forms.)  Anyway, one way to compute this problem, in a specification that encompasses a broad class of functional forms, is to take a finite difference approximation to the partial derivative (see e.g., here).  I do not know if this will give good performance for large $n$, but it is an avenue worth exploring, so I add it here as a suggestion.
Let $\Delta_h$ denote the forward difference operator so that $\Delta_h f(x) = f(x+h) - f(x)$ for a scalar function $f$.  Extend this to a multivariate setting by defining the multivariate operator:
$$\Delta_h^{(k)} \equiv (1, ..., 1, \Delta_h, 1,..., 1),$$
where the finite difference element in the latter vector occurs in the $k$th position, and the unit operators operate on a function multiplying it by one (i.e., leaving it unchanged).  Now, from the first principles definition of Riemann differentiation we have:
$$\frac{\partial}{\partial u_k} C(\mathbf{u}) = \lim_{h \downarrow 0} \frac{\Delta_h^{(k)}}{h} C (\mathbf{u}) \approx \frac{\Delta_h^{(k)}}{h} C (\mathbf{u}) \quad \quad \quad \text{for small } h \approx 0.$$
Hence, if you choose some small value $h \approx 0$, you should get:
$$\begin{equation} \begin{aligned}
\frac{\partial^n C}{\partial u_1 \cdots \partial u_n}(\mathbf{u}) 
= \frac{\partial}{\partial u_1} \cdots \frac{\partial}{\partial u_n} C(\mathbf{u}) 
\approx \frac{\Delta_h^{(1)} \cdots \Delta_h^{(n)}}{h^n} C(\mathbf{u})\\[6pt]
\end{aligned} \end{equation}$$
It should be possible computationally to form the function $\Delta_h^{(1)} \cdots \Delta_h^{(n)} C$ for arbitrary $h$ via combinatorial sums, and since you are using a small values of $h$ you might be able to further approximate this by ignoring some terms.  I'm not sure how computationally-intensive this would be for large $n$, but it is something worth investigating.  That's all I can really add here --- I'm not sure this method would yield a helpful solution, but it is an avenue that I think is worth exploring.
