diff --git a/include/ABACUS_LiebLin.h b/include/ABACUS_LiebLin.h
index 03cda15..12eb717 100644
--- a/include/ABACUS_LiebLin.h
+++ b/include/ABACUS_LiebLin.h
@@ -82,7 +82,7 @@ namespace ABACUS {
     void Set_Free_lambdaocs();
     void Iterate_BAE(DP damping);
     void Iterate_BAE_S(DP damping);
-    void Iterate_BAE_Newton(DP damping);
+    void Iterate_BAE_Newton (DP damping, Vect_DP& RHSBAE, Vect_DP& dlambdaoc, SQMat_DP& Gaudin, Vect_INT& indx);
     void Compute_Energy ();
     void Compute_Momentum ();
     DP Kernel (int a, int b);
diff --git a/src/LIEBLIN/LiebLin_Bethe_State.cc b/src/LIEBLIN/LiebLin_Bethe_State.cc
index 088f559..d2bc934 100644
--- a/src/LIEBLIN/LiebLin_Bethe_State.cc
+++ b/src/LIEBLIN/LiebLin_Bethe_State.cc
@@ -275,8 +275,14 @@ namespace ABACUS {
 
     DP damping = 1.0;
     DP diffsq_prev = 1.0e+6;
+
+    Vect_DP RHSBAE (0.0, N);  // contains RHS of BAEs
+    Vect_DP dlambdaoc (0.0, N);  // contains delta lambdaoc computed from Newton's method
+    SQMat_DP Gaudin (0.0, N);
+    Vect_INT indx (N);
+
     while (diffsq > prec && !is_nan(diffsq) && iter_Newton < 100) {
-      (*this).Iterate_BAE_Newton(damping);
+      (*this).Iterate_BAE_Newton(damping, RHSBAE, dlambdaoc, Gaudin, indx);
       if (diffsq > diffsq_prev && damping > 0.5) damping /= 2.0;
       else if (diffsq < diffsq_prev) damping = 1.0;
       diffsq_prev = diffsq;
@@ -469,12 +475,8 @@ namespace ABACUS {
   /**
      Performs one step of the matrix Newton method on the rapidities.
    */
-  void LiebLin_Bethe_State::Iterate_BAE_Newton (DP damping)
+  void LiebLin_Bethe_State::Iterate_BAE_Newton (DP damping, Vect_DP& RHSBAE, Vect_DP& dlambdaoc, SQMat_DP& Gaudin, Vect_INT& indx)
   {
-    Vect_DP RHSBAE (0.0, N);  // contains RHS of BAEs
-    Vect_DP dlambdaoc (0.0, N);  // contains delta lambdaoc computed from Newton's method
-    SQMat_DP Gaudin (0.0, N);
-    Vect_INT indx (N);
     DP sumtheta = 0.0;
     // for large |lambda|, use atan (lambda) = sgn(lambda) pi/2 - atan(1/lambda)
     int atanintshift = 0;