  In order to accelerate implementation of hyperelastic materials for finite
element analysis, we developed an automatic numerical algorithm that only
requires the strain energy function. This saves the effort on analytical
derivation and coding of stress and tangent modulus, which is time-consuming
and prone to human errors. Using the one-sided Newton difference quotients, the
proposed algorithm first perturbs deformation gradients and calculate the
difference on strain energy to approximate stress. Then, we perturb again to
get difference in stress to approximate tangent modulus. Accuracy of the
approximations were evaluated across the perturbation parameter space, where we
find the optimal amount of perturbation being $10^{-6}$ to obtain stress and
$10^{-4}$ to obtain tangent modulus. Single element verification in ABAQUS with
Neo-Hookean material resulted in a small stress error of only $7\times10^{-5}$
on average across uniaxial compression and tension, biaxial tension and simple
shear situations. A full 3D model with Holzapfel anisotropic material for
artery inflation generated a small relative error of $4\times10^{-6}$ for
inflated radius at $25 kPa$ pressure. Results of the verification tests suggest
that the proposed numerical method has good accuracy and convergence
performance, therefore a good material implementation algorithm in small scale
models and a useful debugging tool for large scale models.
