diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp
index 1f9bc851b92e..98ab059804ce 100644
--- a/Common/include/linear_algebra/CSysSolve.hpp
+++ b/Common/include/linear_algebra/CSysSolve.hpp
@@ -109,8 +109,9 @@ class CSysSolve {
   mutable VectorType p;   /*!< \brief Direction in CG and BCGSTAB. */
   mutable VectorType z;   /*!< \brief Preconditioned residual/direction in CG/BCGSTAB. */
 
-  mutable VectorType r_0; /*!< \brief The "arbitrary" vector in BCGSTAB. */
-  mutable VectorType v;   /*!< \brief BCGSTAB "v" vector (v = A * M^-1 * p). */
+  mutable VectorType r_0;    /*!< \brief The "arbitrary" vector in BCGSTAB. */
+  mutable VectorType v;      /*!< \brief BCGSTAB "v" vector (v = A * M^-1 * p). */
+  mutable VectorType x_best; /*!< \brief BCGSTAB iterate with the smallest residual. */
 
   mutable bool ritz_failed = false;
   mutable unsigned long k = 0, k_new = 0;
diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp
index 8ec1b19b23b3..b480dcacdd95 100644
--- a/Common/src/linear_algebra/CSysSolve.cpp
+++ b/Common/src/linear_algebra/CSysSolve.cpp
@@ -1087,6 +1087,7 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
       p.Initialize(nBlk, nBlkDomain, nVar, nullptr);
       v.Initialize(nBlk, nBlkDomain, nVar, nullptr);
       z.Initialize(nBlk, nBlkDomain, nVar, nullptr);
+      x_best.Initialize(nBlk, nBlkDomain, nVar, nullptr);
 
       bcg_ready = true;
     }
@@ -1135,6 +1136,8 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
   /*--- Initialization ---*/
 
   ScalarType alpha = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0;
+  ScalarType norm_best = norm_r;
+  x_best = x;
   p = ScalarType(0.0);
   v = ScalarType(0.0);
   r_0 = r;
@@ -1195,6 +1198,10 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
       /*--- Check if solution has converged, else output the relative residual if necessary ---*/
 
       norm_r = r.norm();
+      if (norm_r < norm_best) {
+        norm_best = norm_r;
+        x_best = x;
+      }
       if (norm_r < tol * norm0) break;
       if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
         SU2_OMP_MASTER
@@ -1204,6 +1211,13 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
     }
   }
 
+  /*--- BCGSTAB does not reduce the residual monotonically, return the best iterate. ---*/
+
+  if ((config->GetComm_Level() == COMM_FULL) && (norm_r > norm_best)) {
+    x = x_best;
+    norm_r = norm_best;
+  }
+
   /*--- Recalculate final residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {