diff --git a/dumux/material/constraintsolvers/ncpflash.hh b/dumux/material/constraintsolvers/ncpflash.hh
index f1308db95e7bd011335fdb0fdac5c606136d7d11..359350c7f4c490cc065496a5d934d57712b8d207 100644
--- a/dumux/material/constraintsolvers/ncpflash.hh
+++ b/dumux/material/constraintsolvers/ncpflash.hh
@@ -182,10 +182,9 @@ public:
         const int nMax = 50; // <- maximum number of newton iterations
         for (int nIdx = 0; nIdx < nMax; ++nIdx) {
             // regularize fluid state
-            for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
                 if (fluidState.averageMolarMass(phaseIdx) < 1e-10)
                     fluidState.setMoleFraction(phaseIdx, /*compIdx=*/0, 1e-10);
-            }
 
             // calculate Jacobian matrix and right hand side
             linearize_<MaterialLaw>(J, b, fluidState, paramCache, matParams, globalMolarities);
@@ -195,7 +194,23 @@ public:
             // Solve J*x = b
             deltaX = 0;
 
-            try { J.solve(deltaX, b); }
+            try {
+
+                // simple preconditioning of the linear system to reduce condition number
+                int eqIdx = 0;
+                for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                {
+                    for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
+                    {
+                        J[eqIdx] *= 10e-7; // roughly the magnitude of the fugacities
+                        b[eqIdx] *= 10e-7;
+                        ++eqIdx;
+                    }
+                }
+
+                J.solve(deltaX, b);
+
+            }
             catch (Dune::FMatrixError &e)
             {
                 /*
@@ -271,6 +286,11 @@ protected:
             std::cout << fs.density(phaseIdx) << " ";
         std::cout << "\n";
 
+        std::cout << "molar densities: ";
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+            std::cout << fs.molarDensity(phaseIdx) << " ";
+        std::cout << "\n";
+
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
             for (int compIdx = 0; compIdx < numComponents; ++compIdx) {