diff --git a/src/spicelib/devices/bjt/bjtload.c b/src/spicelib/devices/bjt/bjtload.c
index 6e2d832ab..f1ef9439b 100644
--- a/src/spicelib/devices/bjt/bjtload.c
+++ b/src/spicelib/devices/bjt/bjtload.c
@@ -386,6 +386,7 @@ BJTload(GENmodel *inModel, CKTcircuit *ckt)
                         here->BJTtVcrit,&icheck);
                 vbc = DEVpnjlim(vbc,*(ckt->CKTstate0 + here->BJTvbc),vt,
                         here->BJTtVcrit,&ichk1);
+                if (ichk1 == 1) icheck=1;
                 vsub = DEVpnjlim(vsub,*(ckt->CKTstate0 + here->BJTvsub),vt,
                         here->BJTtSubVcrit,&ichk1);
                 if (ichk1 == 1) icheck=1;
@@ -764,7 +765,7 @@ load:
              */
             geqsub = gcsub + gdsub;
             ceqsub=model->BJTtype * model->BJTsubs *
-                    (*(ckt->CKTstate0 + here->BJTcqsub) + cdsub - vsub * geqsub);
+                    (*(ckt->CKTstate0 + here->BJTcqsub) + cdsub - vsub*geqsub);
             ceqbx=model->BJTtype * (*(ckt->CKTstate0 + here->BJTcqbx) -
                     vbx * geqbx);
             ceqbe=model->BJTtype * (cc + cb - vbe * (gm + go + gpi) + vbc *
@@ -779,6 +780,7 @@ load:
                     m * (-ceqbe-ceqbc);
             *(ckt->CKTrhs + here->BJTemitPrimeNode) += m * (ceqbe);
             *(ckt->CKTrhs + here->BJTsubstNode) += m * (-ceqsub);
+
             /*
              *  load y matrix
              */
@@ -786,7 +788,7 @@ load:
             *(here->BJTbaseBasePtr) +=  m * (gx+geqbx);
             *(here->BJTemitEmitPtr) += m * (gepr);
             *(here->BJTcolPrimeColPrimePtr) += m * (gmu+go+gcpr+geqbx);
-            *(here->BJTsubstConSubstConPtr) +=     m * (geqsub);
+            *(here->BJTsubstConSubstConPtr) += m * (geqsub);
             *(here->BJTbasePrimeBasePrimePtr) += m * (gx +gpi+gmu+geqcb);
             *(here->BJTemitPrimeEmitPrimePtr) += m * (gpi+gepr+gm+go);
             *(here->BJTcolColPrimePtr) += m * (-gcpr);
diff --git a/src/spicelib/devices/bjt/bjtpzld.c b/src/spicelib/devices/bjt/bjtpzld.c
index 6f297de82..8226219ac 100644
--- a/src/spicelib/devices/bjt/bjtpzld.c
+++ b/src/spicelib/devices/bjt/bjtpzld.c
@@ -58,8 +58,12 @@ BJTpzLoad(GENmodel *inModel, CKTcircuit *ckt, SPcomplex *s)
             *(here->BJTbaseBasePtr) +=               m * ((gx) + (xcbx) * (s->real));
             *(here->BJTbaseBasePtr + 1) +=           m * ((xcbx) * (s->imag));
             *(here->BJTemitEmitPtr) +=               m * (gepr);
-            *(here->BJTcolPrimeColPrimePtr) +=       m * ((gmu+go+gcpr) + (xcmu+xcsub+xcbx) * (s->real));
-            *(here->BJTcolPrimeColPrimePtr + 1) +=   m * ((xcmu+xcsub+xcbx) * (s->imag));
+            *(here->BJTcolPrimeColPrimePtr) +=       m * ((gmu+go+gcpr) + (xcmu+xcbx) * (s->real));
+            *(here->BJTcolPrimeColPrimePtr + 1) +=   m * ((xcmu+xcbx) * (s->imag));
+
+            *(here->BJTsubstConSubstConPtr) +=       m * (xcsub)* (s->real);
+            *(here->BJTsubstConSubstConPtr + 1) +=   m * (xcsub)* (s->imag);
+
             *(here->BJTbasePrimeBasePrimePtr) +=     m * ((gx+gpi+gmu) + (xcpi+xcmu+xcmcb) * (s->real));
             *(here->BJTbasePrimeBasePrimePtr + 1) += m * ((xcpi+xcmu+xcmcb) * (s->imag));
             *(here->BJTemitPrimeEmitPrimePtr) +=     m * ((gpi+gepr+gm+go) + (xcpi+xgm) * (s->real));
@@ -84,6 +88,10 @@ BJTpzLoad(GENmodel *inModel, CKTcircuit *ckt, SPcomplex *s)
             *(here->BJTemitPrimeBasePrimePtr + 1) += m * ((-xcpi-xgm-xcmcb) * (s->imag));
             *(here->BJTsubstSubstPtr) +=             m * ((xcsub) * (s->real));
             *(here->BJTsubstSubstPtr + 1) +=         m * ((xcsub) * (s->imag));
+            *(here->BJTsubstConSubstPtr) +=          m * ((-xcsub) * (s->real));
+            *(here->BJTsubstConSubstPtr + 1) +=      m * ((-xcsub) * (s->imag));
+            *(here->BJTsubstSubstConPtr) +=          m * ((-xcsub) * (s->real));
+            *(here->BJTsubstSubstConPtr + 1) +=      m * ((-xcsub) * (s->imag));
             *(here->BJTbaseColPrimePtr) +=           m * ((-xcbx) * (s->real));
             *(here->BJTbaseColPrimePtr + 1) +=       m * ((-xcbx) * (s->imag));
             *(here->BJTcolPrimeBasePtr) +=           m * ((-xcbx) * (s->real));