         mat[1][2] = pndcrd[2] - pfccnd[2];
         len = sqrt(mat[1][0]*mat[1][0] + mat[1][1]*mat[1][1]
                  + mat[1][2]*mat[1][2]);
-        mat[1][0] = matinv[0][1] = mat[1][0]/len;
-        mat[1][1] = matinv[1][1] = mat[1][1]/len;
-        mat[1][2] = matinv[2][1] = mat[1][2]/len;
-        mat[2][0] = matinv[0][2] = mat[0][1]*mat[1][2] - mat[0][2]*mat[1][1];
-        mat[2][1] = matinv[1][2] = mat[0][2]*mat[1][0] - mat[0][0]*mat[1][2];
-        mat[2][2] = matinv[2][2] = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
-#else
-        mat[1][0] = matinv[0][1] =  pfcnml[1];
-        mat[1][1] = matinv[1][1] = -pfcnml[0];
+        mat[1][0] = mvt[0][1] = mat[1][0]/len;
+        mat[1][1] = mvt[1][1] = mat[1][1]/len;
+        mat[1][2] = mvt[2][1] = mat[1][2]/len;
+        mat[2][0] = mvt[0][2] = mat[0][1]*mat[1][2] - mat[0][2]*mat[1][1];
+        mat[2][1] = mvt[1][2] = mat[0][2]*mat[1][0] - mat[0][0]*mat[1][2];
+        mat[2][2] = mvt[2][2] = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
+#else
+        mat[1][0] = mvt[0][1] =  pfcnml[1];
+        mat[1][1] = mvt[1][1] = -pfcnml[0];
 #endif
         pdsoln = pidsoln;
-        for (ieq=0; ieq<NEQ; ieq++) {
+        for (ieq=0; ieq<NEQ; ieq+=(NDIM+1)) {
 #if NDIM == 3
             vec[ieq][0] = mat[0][0]*pdsoln[0] + mat[0][1]*pdsoln[1]
                         + mat[0][2]*pdsoln[2];

    int *pfacn, *pfccls;
    double *pidsoln, *pjdsoln, *pdsoln;
    double *pfcnml;
#if NDIM == 3
    int *pfcnds;
    double *pndcrd, *pfccnd;

        pjsoln = exd->soln + jcl*NEQ;
        // rotation and inverse rotation matrices.
        pfcnml = exd->fcnml + ifc*NDIM;
        mat[0][0] = mvt[0][0] = pfcnml[0];
        mat[0][1] = mvt[1][0] = pfcnml[1];
#if NDIM == 3
        mat[0][2] = mvt[2][0] = pfcnml[2];
        pfcnds = exd->fcnds + ifc*(FCMND+1);
        pndcrd = exd->ndcrd + pfcnds[1]*NDIM;
        pfccnd = exd->fccnd + ifc*NDIM;
