Best strategy for verification of tangent for hyperelastic material

[WayneKagawa]

Let's say I have a strain energy density function W = W(I1,I2) + W(I3).
The second PK stress is then defined as S_{ij} = ∂W/∂E_{ij} and the constitutive tensor C_{ijkl} = ∂/∂E_{kl}(∂W/∂E_{ij}).

If I want to verify that I'm deriving the correct analytical constitutive and stress tensor equations, what's the best way without using a commercial code?

The strategy I chose was using a MATLAB code to simulate uni-axial tension the following way:
lambda = [1.0 -> 1.5] (100x1 array of applied stretch to material)
For i in lambda:

Define deformation gradient F​

Get right cauchy green tensor C and Lagrange strain E​

Get invariants I1, I2 and I3 of C​

Get S^(analytical)_{ij}=∂W/∂E_{ij} by plugging in invariants​

Get stress from product of constitutive tensor and lagrange strain, S^(tangent)_{ij} = C_{ijkl}E_{kl}​

Then compare the two stresses, S^(analytical)_{ij} and S^(tangent)_{ij}(either through a plot or through a point wise error).

Is this way correct? Essentially, is S_{ij} = C_{ijkl}E_{kl} valid at all points in the range of stretches for hyper-elastic materials?

Thank you!