diff --git a/src/TBA/TBA_LiebLin.cc b/src/TBA/TBA_LiebLin.cc
index 0a3c4b2..a25947a 100644
--- a/src/TBA/TBA_LiebLin.cc
+++ b/src/TBA/TBA_LiebLin.cc
@@ -643,37 +643,55 @@ namespace ABACUS {
     return(sbar);
   }
 
+
   LiebLin_Bethe_State Discretized_LiebLin_Bethe_State (DP c_int, DP L, int N, const Root_Density& rho)
   {
     // This function returns the Bethe state at finite size which is
     // the closest approximation to the continuum density rho(lambda)
 
-    // Each time N \int_{-\infty}^\lambda d\lambda' \rho(\lambda') crosses a half integer, add a particle:
+    // Check that the provided rho has the expected filling
+    DP rho_check = 0.0;
+    for (int i = 0; i < rho.Npts; ++i) rho_check += L * rho.value[i] * rho.dlambda[i];
+    // The integral of rho should be between N - 0.5 and N + 0.5 for the algorithm to work
+    if (fabs(rho_check - N) > 0.5)
+      cout << "WARNING: integral of rho != N in Discretized_LiebLin_Bethe_State, "
+	   << "consistent discretization is impossible. \int rho dlambda = "
+	   << rho_check << ", N = " << N << ", L = " << L << ", N/L = " << N/L << endl;
+
+    // Using the counting function
+    // c(\lambda) \equiv L \int_{-\infty}^\lambda d\lambda' \rho(\lambda'),
+    // we seek to find lambda[i] such that c(lambda[i]) = i + 0.5, i = 0...N-1
     DP integral = 0.0;
     DP integral_prev = 0.0;
     int Nfound = 0;
-    Vect<DP> lambda_found(0.0, 2*N);
+    Vect<DP> lambda_found(0.0, N);
+    int nr_to_set = 0;
 
     for (int i = 0; i < rho.Npts; ++i) {
+      // Note: integral_prev is the counting function
+      // at rapidity rho.lambda[i-1] + 0.5* rho.dlambda[i-1]
+      // which equals rho.lambda[i] - 0.5* rho.dlambda[i]
       integral_prev = integral;
+      // Note: integral gives the value of the counting function
+      // at rapidity rho.lambda[i] + 0.5* rho.dlambda[i]
       integral += L * rho.value[i] * rho.dlambda[i];
-      if (integral > Nfound + 0.5) {
-	// Subtle error: if the rho is too discontinuous, i.e. if more than one rapidity is found, must correct for this.
-	if (integral > Nfound + 1.5 && integral < Nfound + 2.5) { // found two rapidities
-	  lambda_found[Nfound++] = 0.25 * (3.0 * rho.lambda[i-1] + rho.lambda[i]);
-	  lambda_found[Nfound++] = 0.25 * (rho.lambda[i-1] + 3.0 * rho.lambda[i]);
-	}
-	else {
-	  // Better: center the lambda_found between these points:
-	  lambda_found[Nfound] = 0.5 * (rho.lambda[i-1] + rho.lambda[i]);
-	  Nfound++;
-	}
+      // We already have filled the RHS of the counting equation up to Nfound - 0.5.
+      // Therefore, the additional number found is
+      nr_to_set = floor(integral + 0.5 - Nfound);
+      // if (nr_to_set > 1)
+      // 	cout << "WARNING: setting " << nr_to_set << " rapidities in one step in Discretized_LiebLin_Bethe_State" << endl;
+      for (int n = 1; n <= nr_to_set; ++n) {
+	// Solve c(lambda[Nfound]) = Nfound + 0.5 for lambda[Nfound],
+	// namely (using linear interpolation)
+	// integral_prev + (lambda - rho.lambda[i-1] - 0.5*rho.dlambda[i-1]) * (integral - integra_prev)/rho.dlambda[i] = Nfound + 0.5
+	lambda_found[Nfound] = rho.lambda[i-1] + 0.5*rho.dlambda[i-1] + (Nfound + 0.5 - integral_prev) * rho.dlambda[i]/(integral - integral_prev);
+	Nfound++;
       }
     }
 
     Vect<DP> lambda(N);
     // Fill up the found rapidities:
-    for (int il = 0; il < ABACUS::min(N, Nfound); ++il) lambda[il] = lambda_found[il];
+    for (int il = 0; il < abacus::min(N, Nfound); ++il) lambda[il] = lambda_found[il];
     // If there are missing ones, put them at the end; ideally, this should never be called
     for (int il = Nfound; il < N; ++il) lambda[il] = lambda_found[Nfound-1] + (il - Nfound + 1) * (lambda_found[Nfound-1] - lambda_found[Nfound-2]);
 
@@ -687,6 +705,9 @@ namespace ABACUS {
       Ix2[i] = 2.0 * floor((L* lambda[i] + sum)/twoPI + 0.5 * (N%2 ? 1 : 2)) + (N%2) - 1;
     }
 
+    // Check that the Ix2 are all ordered
+    for (int i = 0; i < N-1; ++i) if (Ix2[i] >= Ix2[i+1]) cout << "Alert: Ix2 not ordered around index i = " << i << ": Ix2[i] = " << Ix2[i] << "\tIx2[i+1] = " << Ix2[i+1] << endl;
+
     // Check that the quantum numbers are all distinct:
     bool allOK = false;
     while (!allOK) {