Plane Stress (Shell) incompressible hyperelasticity VUMAT written in C -- Example & Instability problem

Based on examples from Dr. Bergstorm, I wrote a VUMAT in C to implement Yeoh hyperelasticity for both 3D and plane stress whereas plane stress assumes incompressible material.  The incompressibility is enforced by F33 = 1./(F11*F22-F12*F21) and plane stress by finding duDi1 so that S33=0. 

The 3D version works well and results match ABQ almost exactly including time increment selected by ABQ for both built-in material model and vumat. 

However, the plane stress version still has problem as below, and was hoping getting some help here. 

When using small fixed time step or use "element by element" option for increment size control, vumat matches ABQ results almost exactly. in this case, vumat is compared with ABQ parallel network model with only the elastic network active as it also assumes incompressible in plane stress.  However, if use automatic time increment control without "element by element" option, vumat goes into instable fairly quickly and almost all of sudden in a few increments after performing well for thousands of increments. I just could not figure out what the problem is in the vumat.    it seems that vumat drives ABQ to choose too large increment resulting in instability. 

files can be downloaded here

https://1drv.ms/u/s!AvzF3H-Ljv02g8FsQCGROQuwyuleew?e=E2WWmP

elastic-D.cpp -- vumat  written in c. should work with abaqus 2019 - 2021, double precision

shell-shear-ABQ.inp  pure shear case, ABQ built-in material model

shell-shear-vumat.inp pure shear case, vumat

shell-uniaxial-ABQ.inp uniaxial tensile case, ABQ built-in material model

shell-uniaxial-vumat.inp uniaxial tensile case, vumat