diff --git a/apps/regression/model/model.cpp b/apps/regression/model/model.cpp
index def197f1a..42dc50a2e 100644
--- a/apps/regression/model/model.cpp
+++ b/apps/regression/model/model.cpp
@@ -65,7 +65,9 @@ void Model::fitLevenbergMarquardt(Store * store, int series, double * modelCoeff
       operandsB[j] = betaCoefficient(store, series, modelCoefficients, j);
     }
     double modelCoefficientSteps[Model::k_maxNumberOfCoefficients];
-    solveLinearSystem(modelCoefficientSteps, coefficientsAPrime, operandsB, n, context);
+    if (solveLinearSystem(modelCoefficientSteps, coefficientsAPrime, operandsB, n, context) < 0) {
+      break;
+    }
     // Compare new chi2 with the previous value
     double newModelCoefficients[Model::k_maxNumberOfCoefficients];
     for (int i = 0; i < n; i++) {
@@ -146,7 +148,7 @@ double Model::betaCoefficient(Store * store, int series, double * modelCoefficie
   return result;
 }
 
-void Model::solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Context * context) {
+int Model::solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Context * context) {
   int n = solutionDimension;
   assert(n <= k_maxNumberOfCoefficients);
   double coefficientsSave[k_maxNumberOfCoefficients * k_maxNumberOfCoefficients];
@@ -154,18 +156,25 @@ void Model::solveLinearSystem(double * solutions, double * coefficients, double
     coefficientsSave[i] = coefficients[i];
   }
   int inverseResult = Matrix::ArrayInverse(coefficients, n, n);
-  if (inverseResult < 0) {
+  int numberOfMatrixModifications = 0;
+  while (inverseResult < 0 && numberOfMatrixModifications < k_maxMatrixInversionFixIterations) {
     // If the matrix was not invertible, modify it a bit
     for (int i = 0; i < n; i ++) {
       coefficientsSave[i*n+i] = (1 + ((double)i)/((double)n)) * coefficientsSave[i*n+i];
     }
-    int inverseSecondResult = Matrix::ArrayInverse(coefficientsSave, n, n);
-    assert(inverseSecondResult >= 0);
+    inverseResult = Matrix::ArrayInverse(coefficientsSave, n, n);
+    numberOfMatrixModifications++;
+  }
+  if (inverseResult < 0) {
+    return - 1;
+  }
+  if (numberOfMatrixModifications > 0) {
     for (int i = 0; i < n*n; i++) {
       coefficients[i] = coefficientsSave[i];
     }
   }
   Multiplication::computeOnArrays<double>(coefficients, constants, solutions, n, n, 1);
+  return 0;
 }
 
 }
diff --git a/apps/regression/model/model.h b/apps/regression/model/model.h
index 71e788d05..fd91b4b5b 100644
--- a/apps/regression/model/model.h
+++ b/apps/regression/model/model.h
@@ -42,6 +42,7 @@ private:
 
   // Levenberg-Marquardt
   static constexpr double k_maxIterations = 1000;
+  static constexpr double k_maxMatrixInversionFixIterations = 10;
   static constexpr double k_initialLambda = 0.001;
   static constexpr double k_lambdaFactor = 10;
   static constexpr double k_chi2ChangeCondition = 0.001;
@@ -52,7 +53,7 @@ private:
   double alphaPrimeCoefficient(Store * store, int series, double * modelCoefficients, int k, int l, double lambda) const;
   double alphaCoefficient(Store * store, int series, double * modelCoefficients, int k, int l) const;
   double betaCoefficient(Store * store, int series, double * modelCoefficients, int k) const;
-  void solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Poincare::Context * context);
+  int solveLinearSystem(double * solutions, double * coefficients, double * constants, int solutionDimension, Poincare::Context * context);
 };
 
 }