This is an attempt which I didn't completely work through, but too long for the comments section. It might be useful to put it here as another basic alternative for very low $k$. It does not require explicit differentiation + MCMC (but it does perform numerical differentiation, without MCMC).
Algorithm
For small $\varepsilon > 0$:
- Draw $u_1 \sim C_1 \equiv C(U_1 = u_1,U_2 = 1,\ldots, U_k = 1)$. This can be easily done by drawing $\eta \sim \text{Uniform}[0,1]$ and computing $C_1^{-1}(\eta)$ (which, if anything, can be easily done numerically). This is a draw from the marginal pdf $u_1 \sim \kappa(u_1)$.
- For $j = 2\ldots k$
- Define $$D_j^{(\varepsilon)}(u_j|u_1,\ldots,u_{j-1}) \equiv \Pr\left( u_1 - \frac{\varepsilon}{2} \le U_1 \le u_1 + \frac{\varepsilon}{2} \land \dots \land u_{j-1} - \frac{\varepsilon}{2} \le U_{j-1} \le u_{j-1} + \frac{\varepsilon}{2} \land U_{j} \le u_j \land U_{j+1} \le 1 \dots \land U_{k} \le 1\right),$$ which can be computed as a difference of $C$ evaluated at various points (which in the naive way needs $O(2^{j-1})$ evaluations of $C$ for every evaluation of $D_j^{(\varepsilon)}$). $D_j^{(\varepsilon)}$ is the $\varepsilon$-approximate marginal conditional of $u_j$ given $u_1, \ldots, u_{j-1}$.
- Draw $u_j \sim D_j^{(\varepsilon)}(u_j|u_1,\ldots,u_{j-1})$ as per point 1, which again should be easy to do with numerical inversion.
Discussion
This algorithm should generate i.i.d. samples from an $\varepsilon$-approximation of $C(u_1,\ldots,u_k)$, where $\varepsilon$ merely depends on numerical precision. There are practical technicalities to refine the approximation and make it numerically stable.
The obvious problem is that computational complexity scales as $O(2^{k})$, so, to put it generously, this is not very general in terms of $k$ (but the example you linked had $k = 3$, so perhaps this method is not completely useless -- I am not familiar with the typical scenario in which you would have access to the cdf). On the other hand, for very low-dimensional distributions it could work, and the cost is compensated by the fact that, unlike the other generic solution of "differentiating + MCMC", there is no need to compute derivatives, samples are i.i.d. and there is no tuning (aside the choice of $\varepsilon$, which should just be something slightly above machine precision). And perhaps there are ways to make this better than the naive approach.
As I mentioned, this is off the top of my head so there might be other issues.