diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index 5911d86d87ae92b04137d5566dc146559ea2b76a..253e6369b9f0d4d087e7cc8eab14c70b3b04230e 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -379,29 +379,70 @@ class PolynomialP
       }
       return os;
     } else if constexpr (Dimension == 2) {
-      for (size_t j = N; j > 0; --j) {
-        for (size_t i = N - j; i > 0; --i) {
-          double coef = P[i, j];
-          if (coef != 0.) {
-            if (coef > 0.) {
-              os << " + ";
-            } else if (coef < 0.) {
-              os << " - ";
-            }
-            if (coef != 1 && coef != -1) {
-              os << std::abs(coef);
+      size_t i = 0;
+      size_t j = N;
+      TinyVector<Dimension, size_t> rel_pos(i, j);
+      double coef = P[rel_pos];
+      if (coef != 0.) {
+        if (coef < 0.) {
+          os << " - ";
+        }
+        if (coef != 1 && coef != -1) {
+          os << std::abs(coef);
+        }
+        os << "y^" << N;
+      }
+      size_t degree = N;
+      for (size_t k = size_coef - 1; k > 0; --k) {
+        if (j > 0) {
+          j--;
+          i++;
+        } else {
+          degree--;
+          j = degree;
+          i = 0;
+        }
+        rel_pos     = TinyVector<Dimension, size_t>(i, j);
+        double coef = P[rel_pos];
+        if (coef != 0.) {
+          if (coef > 0.) {
+            os << " + ";
+          } else if (coef < 0.) {
+            os << " - ";
+          }
+          if ((coef != 1 && coef != -1) || (i == 0 && j == 0)) {
+            os << std::abs(coef);
+          }
+          if (i == 0 && j == 0)
+            continue;
+          if (i == 0) {
+            if (j != 1) {
+              os << "y^" << j;
+            } else {
+              os << "y";
             }
-            if (i == 0) {
-              os << "y" << j;
-            } else if (j == 0) {
+          } else if (j == 0) {
+            if (i == 1) {
+              os << "x";
+            } else {
               os << "x^" << i;
+            }
+          } else {
+            if (i == 1 && j == 1) {
+              os << "xy";
+            } else if (i == 1) {
+              os << "x"
+                 << "y^" << j;
+            } else if (j == 1) {
+              os << "x^" << i << "y";
             } else {
               os << "x^" << i << "y^" << j;
             }
-            all_coef_zero = false;
           }
+          all_coef_zero = false;
         }
       }
+
       return os;
     } else {
       throw NotImplementedError("Not yet Available in 3D");