In this study, we presented a numerically robust procedure to evaluate 4th order tangent moduli which are vital for acheiving quadratic convergence of global Newton-Raphson scheme.
In this study, the proposed method is verified for hyperelastic models alone. However, the same can be extended to other constitutive models.
The paper can be found at
Numerically approximated Cauchy integral (NACI) for implementation of constitutive models
10.1016/j.finel.2014.05.016
http://www.sciencedirect.com/science/article/pii/S0168874X14001061
A side note
The procedure provided in this paper can also be used to evaluate numerically accurate (see paper for definition) very high order derivatives of differentiable functions.