gh1h3h3h3=-3*(L4*R12*R31^2*R32 - cb*IL5*R13*R31^2*R32 + L4*R11*R31*R32^2 + 
  L2*R12*R32^3 - cb*IL5*R12*R31^2*R33 + cb^2*L4*R13*R31^2*R33 - 
  2*cb*IL5*R11*R31*R32*R33 + cb^2*L2*R13*R32^2*R33 + cb^2*L4*R11*R31*R33^2 + 
  cb^2*L2*R12*R32*R33^2 + cb^4*L2*R13*R33^3 + R12*R31^2*R32*RL5 + 
  R11*R31*R32^2*RL5 - cb^2*R13*R31^2*R33*RL5 - cb^2*R11*R31*R33^2*RL5 - 
  IL5*R13*R31*R32^2*sb - 2*IL5*R12*R31*R32*R33*sb - IL5*R11*R32^2*R33*sb + 
  3*cb^2*IL5*R13*R31*R33^2*sb + cb^2*IL5*R11*R33^3*sb - 
  4*cb*R13*R31*R32*R33*RL5*sb - 2*cb*R12*R31*R33^2*RL5*sb - 
  2*cb*R11*R32*R33^2*RL5*sb + L4*R13*R32^2*R33*sb^2 + L4*R12*R32*R33^2*sb^2 + 
  3*cb*IL5*R13*R32*R33^2*sb^2 + cb*IL5*R12*R33^3*sb^2 + 
  2*cb^2*L4*R13*R33^3*sb^2 - R13*R32^2*R33*RL5*sb^2 - 
  R12*R32*R33^2*RL5*sb^2 + 2*cb^2*R13*R33^3*RL5*sb^2 + 
  L1*(R11*R31 + R13*R33*sb^2)*(R31^2 + R33^2*sb^2) + 
  L3*(R11*R31*(R32^2 + cb^2*R33^2) + R12*R32*(R31^2 + R33^2*sb^2) + 
    R13*R33*(R32^2*sb^2 + cb^2*(R31^2 + 2*R33^2*sb^2))))
