=nothing major

src/residue.cpp

@@ -1789,6 +1789,7 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data)
 
     /*implement res @ left boundary i.e., first grid point*/
     Rres(1) = Rdot(1);
+    sum = ZERO ;
     if (dropType == 0) {
         // Yres(1,data->k_drop[0]) = Ydot(1,data->k_drop[0]);
         for (size_t k = 1; k <= nsp; k++) {
@@ -1798,12 +1799,18 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data)
         for (size_t k = 1; k <= nsp; k++) {
             if (k == k_drop[0]) {
                 Yres(1, k) = Y(2, k) - Y(1, k);
-            } else if (k == k_drop[1]) {
-                Yres(1, k) = ONE - Y(1, k_drop[0]) - Y(1, k);
-            } else {
+                sum += Y(1, k);
+            }
+//            else if (k == k_drop[1]) {
+//                //Yres(1, k) = ONE - Y(1, k_drop[0]) - Y(1, k);
+//                Yres(1,k) =
+//            }
+            else if (k != k_drop[1] && k != k_drop[0]) {
                 Yres(1, k) = Ydot(1, k);
+                sum += Y(1, k);
             }
         }
+        Yres(1, k_drop[1]) = ONE - sum - Y(1, k_drop[1]);
     }
     Pres(1) = P(2) - P(1);
     Mdotres(1) = Mdot(2) - Mdot(1);