diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index 42f0c46663019d7f4e65220c8e25fbd6a89c94e6..76f6b46d63512f5bb32f2a6e0dc466a48e03e627 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <analysis/CubeGaussQuadrature.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
 #include <analysis/GaussQuadratureDescriptor.hpp>
 #include <analysis/QuadratureManager.hpp>
 #include <analysis/SquareGaussQuadrature.hpp>
@@ -616,7 +617,8 @@ integrate(const PolynomialP<N, Dimension> Q, const std::array<TinyVector<Dimensi
   } else if constexpr (Dimension == 2) {
     static_assert(P <= 4, "In 2D number of positions should be lesser or equal to 4");
     if constexpr (P == 2) {
-      const QuadratureFormula<1>& lN = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(N));
+      const QuadratureFormula<1>& lN =
+        QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(N + 1));
       const LineTransformation<2> t{positions[0], positions[1]};
       double value = 0.;
       for (size_t i = 0; i < lN.numberOfPoints(); ++i) {
@@ -635,9 +637,10 @@ integrate(const PolynomialP<N, Dimension> Q, const std::array<TinyVector<Dimensi
       }
       integral = t.jacobianDeterminant() * value;
     } else {
-      const QuadratureFormula<2>& lN = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(N));
-      auto point_list                = lN.pointList();
-      auto weight_list               = lN.weightList();
+      const QuadratureFormula<2>& lN =
+        QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(N + 1));
+      auto point_list  = lN.pointList();
+      auto weight_list = lN.weightList();
       SquareTransformation<2> s{positions[0], positions[1], positions[2], positions[3]};
       auto value = weight_list[0] * s.jacobianDeterminant(point_list[0]) * Q(s(point_list[0]));
       for (size_t i = 1; i < weight_list.size(); ++i) {
diff --git a/src/analysis/TaylorPolynomial.hpp b/src/analysis/TaylorPolynomial.hpp
index 1406a5989a4393e90933f4e59ff232f793311082..19ceec686a5f18728506f3daba9e27daa888fbcb 100644
--- a/src/analysis/TaylorPolynomial.hpp
+++ b/src/analysis/TaylorPolynomial.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <analysis/CubeGaussQuadrature.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
 #include <analysis/GaussQuadratureDescriptor.hpp>
 #include <analysis/QuadratureManager.hpp>
 #include <analysis/SquareGaussQuadrature.hpp>
@@ -169,7 +170,7 @@ class TaylorPolynomial
 
   PUGS_INLINE
   constexpr double
-  operator[](const TinyVector<Dimension, size_t> relative_pos) const
+  operator[](const TinyVector<Dimension, size_t>& relative_pos) const
   {
     size_t total_degree = 0;
     for (size_t i = 0; i < Dimension; ++i) {
@@ -179,8 +180,7 @@ class TaylorPolynomial
     }
     Assert((total_degree <= N), "The sum of the degrees of the coefficient you are looking for is greater than the "
                                 "degree of the polynomial");
-    TinyVector<size_coef> absolute_coefs = this->coefficients();
-    size_t absolute_position             = 0;
+    size_t absolute_position = 0;
     if constexpr (Dimension == 1) {
       absolute_position = relative_pos[0];
     } else if constexpr (Dimension == 2) {
@@ -195,12 +195,12 @@ class TaylorPolynomial
                           total_sub_degree * (total_sub_degree + 1) / 2 + relative_pos[2];
     }
 
-    return absolute_coefs[absolute_position];
+    return m_coefficients[absolute_position];
   }
 
   PUGS_INLINE
-  constexpr double
-  operator[](const TinyVector<Dimension, size_t> relative_pos)
+  constexpr double&
+  operator[](const TinyVector<Dimension, size_t>& relative_pos)
   {
     size_t total_degree = 0;
     for (size_t i = 0; i < Dimension; ++i) {
@@ -210,8 +210,7 @@ class TaylorPolynomial
     }
     Assert((total_degree <= N), "The sum of the degrees of the coefficient you are looking for is greater than the "
                                 "degree of the polynomial");
-    TinyVector<size_coef> absolute_coefs = this->coefficients();
-    size_t absolute_position             = 0;
+    size_t absolute_position = 0;
     if constexpr (Dimension == 1) {
       absolute_position = relative_pos[0];
     } else if constexpr (Dimension == 2) {
@@ -226,7 +225,7 @@ class TaylorPolynomial
                           total_sub_degree * (total_sub_degree + 1) / 2 + relative_pos[2];
     }
 
-    return absolute_coefs[absolute_position];
+    return m_coefficients[absolute_position];
   }
 
   PUGS_INLINE
@@ -525,12 +524,10 @@ class TaylorPolynomial
                  << "(y - " << P.x0()[1] << ")^" << j;
             } else if (j == 1) {
               os << "(x - " << P.x0()[0] << ")"
-                 << "^i"
-                 << "(y - " << P.x0()[1] << ")";
+                 << "^" << i << "(y - " << P.x0()[1] << ")";
             } else {
               os << "(x - " << P.x0()[0] << ")"
-                 << "^i"
-                 << "(y - " << P.x0()[1] << ")^" << j;
+                 << "^" << i << "(y - " << P.x0()[1] << ")^" << j;
             }
           }
           all_coef_zero = false;
@@ -663,9 +660,10 @@ integrate(const TaylorPolynomial<N, Dimension> Q, const std::array<TinyVector<Di
       }
       integral = t.jacobianDeterminant() * value;
     } else {
-      const QuadratureFormula<2>& lN = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(N));
-      auto point_list                = lN.pointList();
-      auto weight_list               = lN.weightList();
+      const QuadratureFormula<2>& lN =
+        QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(N + 1));
+      auto point_list  = lN.pointList();
+      auto weight_list = lN.weightList();
       SquareTransformation<2> s{positions[0], positions[1], positions[2], positions[3]};
       auto value = weight_list[0] * s.jacobianDeterminant(point_list[0]) * Q(s(point_list[0]));
       for (size_t i = 1; i < weight_list.size(); ++i) {
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
index 7b46193c38434c2400a34fa0248f143461a2adac..bd6d4a82d57f16447ffdb70bba302eee1e9267e5 100644
--- a/tests/test_PolynomialP.cpp
+++ b/tests/test_PolynomialP.cpp
@@ -1,5 +1,6 @@
 #include <Kokkos_Core.hpp>
 #include <algebra/TinyMatrix.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
 #include <analysis/GaussQuadratureDescriptor.hpp>
 #include <analysis/PolynomialP.hpp>
 #include <analysis/QuadratureManager.hpp>
@@ -167,10 +168,12 @@ TEST_CASE("PolynomialP", "[analysis]")
   SECTION("integrate")
   {
     TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+    TinyVector<10> coef2(1, -2, 10, 7, 2, 9, 0, 0, 0, 1);
     std::array<TinyVector<2>, 4> positions;
     std::array<TinyVector<2>, 3> positions2;
     std::array<TinyVector<2>, 2> positions4;
     std::array<TinyVector<2>, 4> positions3;
+    std::array<TinyVector<2>, 4> positions5;
     positions[0]  = TinyVector<2>{0, 0};
     positions[1]  = TinyVector<2>{0, 0.5};
     positions[2]  = TinyVector<2>{0.3, 0.7};
@@ -184,21 +187,40 @@ TEST_CASE("PolynomialP", "[analysis]")
     positions3[1] = TinyVector<2>{1, 0};
     positions3[2] = TinyVector<2>{1, 1};
     positions3[3] = TinyVector<2>{0, 1};
+    positions5[0] = TinyVector<2>{0, 0};
+    positions5[1] = TinyVector<2>{1, 0};
+    positions5[2] = TinyVector<2>{1, 3};
+    positions5[3] = TinyVector<2>{1, 1};
 
     PolynomialP<2, 2> P(coef);
+    PolynomialP<3, 2> Q(coef2);
     auto p1 = [](const TinyVector<2>& X) {
       const double x = X[0];
       const double y = X[1];
       return 1 - 2. * x + 10 * y + 7 * x * x + 2 * x * y + 9 * y * y;
     };
-    const QuadratureFormula<2>& l3 = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(3));
-    auto point_list                = l3.pointList();
-    auto weight_list               = l3.weightList();
+    auto q1 = [](const TinyVector<2>& X) {
+      const double x = X[0];
+      const double y = X[1];
+      return 1 - 2. * x + 10 * y + 7 * x * x + 2 * x * y + 9 * y * y + y * y * y;
+    };
+    const QuadratureFormula<2>& l3 =
+      QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(4));
+    auto point_list  = l3.pointList();
+    auto weight_list = l3.weightList();
     SquareTransformation<2> s{positions[0], positions[1], positions[2], positions[3]};
-    auto value = weight_list[0] * s.jacobianDeterminant(point_list[0]) * p1(s(point_list[0]));
+    auto value = weight_list[0] * s.jacobianDeterminant(point_list[0]) * q1(s(point_list[0]));
     for (size_t i = 1; i < weight_list.size(); ++i) {
-      value += weight_list[i] * s.jacobianDeterminant(point_list[i]) * p1(s(point_list[i]));
+      value += weight_list[i] * s.jacobianDeterminant(point_list[i]) * q1(s(point_list[i]));
     }
+    // const QuadratureFormula<2>& l3bis =
+    //   QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(4));
+    // auto point_listbis  = l3bis.pointList();
+    // auto weight_listbis = l3bis.weightList();
+    // auto valuebis       = weight_listbis[0] * s.jacobianDeterminant(point_listbis[0]) * p1(s(point_listbis[0]));
+    // for (size_t i = 1; i < weight_listbis.size(); ++i) {
+    //   valuebis += weight_listbis[i] * s.jacobianDeterminant(point_listbis[i]) * p1(s(point_listbis[i]));
+    // }
     const QuadratureFormula<2>& t3 = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(3));
     auto point_list2               = t3.pointList();
     auto weight_list2              = t3.weightList();
@@ -208,7 +230,7 @@ TEST_CASE("PolynomialP", "[analysis]")
       value2 += weight_list2[i] * p1(t(point_list2[i]));
     }
     value2 *= t.jacobianDeterminant();
-    const QuadratureFormula<1>& l1 = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+    const QuadratureFormula<1>& l1 = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(2));
     const LineTransformation<2> u{positions4[0], positions4[1]};
     double value4 = 0.;
     for (size_t i = 0; i < l1.numberOfPoints(); ++i) {
@@ -220,10 +242,12 @@ TEST_CASE("PolynomialP", "[analysis]")
       value3 += weight_list[i] * s3.jacobianDeterminant(point_list[i]) * p1(s3(point_list[i]));
     }
 
-    REQUIRE(value == Catch::Approx(integrate(P, positions)));
+    REQUIRE(value == Catch::Approx(integrate(Q, positions)));
     REQUIRE(value2 == Catch::Approx(integrate(P, positions2)));
     REQUIRE(value3 == Catch::Approx(integrate(P, positions3)));
     REQUIRE(value4 == Catch::Approx(integrate(P, positions4)));
+    // REQUIRE(valuebis == Catch::Approx(integrate(P, positions)));
+    //    REQUIRE(valuebis == value);
   }
 
   // // SECTION("product")
diff --git a/tests/test_TaylorPolynomial.cpp b/tests/test_TaylorPolynomial.cpp
index b5871d246d4841c50e55f4618f93d731c7fe1ea4..4d9149232202adf49de64107fb3bfbe8e20eed22 100644
--- a/tests/test_TaylorPolynomial.cpp
+++ b/tests/test_TaylorPolynomial.cpp
@@ -160,6 +160,7 @@ TEST_CASE("TaylorPolynomial", "[analysis]")
   SECTION("integrate")
   {
     TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+    TinyVector<10> coef2(1, -2, 10, 7, 2, 9, 0, 0, 0, 1);
     std::array<TinyVector<2>, 4> positions;
     std::array<TinyVector<2>, 3> positions2;
     std::array<TinyVector<2>, 4> positions3;
@@ -178,41 +179,60 @@ TEST_CASE("TaylorPolynomial", "[analysis]")
     positions4[0] = TinyVector<2>{0, 0.5};
     positions4[1] = TinyVector<2>{0.3, 0.7};
     TaylorPolynomial<2, 2> P(coef, x0);
+    TaylorPolynomial<2, 2> Q(coef, zero);
+    TaylorPolynomial<3, 2> R(coef2, x0);
     auto p1 = [](const TinyVector<2>& X) {
+      const double x = X[0];
+      const double y = X[1];
+      return 1 - 2. * (x - 1.) + 10 * (y + 1) + 7 * (x - 1.) * (x - 1) + 2 * (x - 1) * (y + 1) + 9 * (y + 1) * (y + 1) +
+             (y + 1) * (y + 1) * (y + 1);
+    };
+    auto p2 = [](const TinyVector<2>& X) {
       const double x = X[0];
       const double y = X[1];
       return 1 - 2. * (x - 1.) + 10 * (y + 1) + 7 * (x - 1.) * (x - 1) + 2 * (x - 1) * (y + 1) + 9 * (y + 1) * (y + 1);
     };
-    const QuadratureFormula<2>& l3 = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(3));
-    auto point_list                = l3.pointList();
-    auto weight_list               = l3.weightList();
+    auto q1 = [](const TinyVector<2>& X) {
+      const double x = X[0];
+      const double y = X[1];
+      return 1 - 2. * x + 10 * y + 7 * x * x + 2 * x * y + 9 * y * y;
+    };
+    const QuadratureFormula<2>& l3 =
+      QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(4));
+    auto point_list  = l3.pointList();
+    auto weight_list = l3.weightList();
     SquareTransformation<2> s{positions[0], positions[1], positions[2], positions[3]};
     auto value = weight_list[0] * s.jacobianDeterminant(point_list[0]) * p1(s(point_list[0]));
     for (size_t i = 1; i < weight_list.size(); ++i) {
       value += weight_list[i] * s.jacobianDeterminant(point_list[i]) * p1(s(point_list[i]));
     }
+    auto valuebis = weight_list[0] * s.jacobianDeterminant(point_list[0]) * q1(s(point_list[0]));
+    for (size_t i = 1; i < weight_list.size(); ++i) {
+      valuebis += weight_list[i] * s.jacobianDeterminant(point_list[i]) * q1(s(point_list[i]));
+    }
     SquareTransformation<2> s3{positions3[0], positions3[1], positions3[2], positions3[3]};
-    auto value3 = weight_list[0] * s3.jacobianDeterminant(point_list[0]) * p1(s3(point_list[0]));
+    auto value3 = weight_list[0] * s3.jacobianDeterminant(point_list[0]) * p2(s3(point_list[0]));
     for (size_t i = 1; i < weight_list.size(); ++i) {
-      value3 += weight_list[i] * s3.jacobianDeterminant(point_list[i]) * p1(s3(point_list[i]));
+      value3 += weight_list[i] * s3.jacobianDeterminant(point_list[i]) * p2(s3(point_list[i]));
     }
     const QuadratureFormula<2>& t3 = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(3));
     auto point_list2               = t3.pointList();
     auto weight_list2              = t3.weightList();
     TriangleTransformation<2> t{positions2[0], positions2[1], positions2[2]};
-    auto value2 = weight_list2[0] * p1(t(point_list2[0]));
+    auto value2 = weight_list2[0] * p2(t(point_list2[0]));
     for (size_t i = 1; i < weight_list2.size(); ++i) {
-      value2 += weight_list2[i] * p1(t(point_list2[i]));
+      value2 += weight_list2[i] * p2(t(point_list2[i]));
     }
     value2 *= t.jacobianDeterminant();
     const QuadratureFormula<1>& l1 = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
     const LineTransformation<2> u{positions4[0], positions4[1]};
     double value4 = 0.;
     for (size_t i = 0; i < l1.numberOfPoints(); ++i) {
-      value4 += l1.weight(i) * u.velocityNorm() * p1(u(l1.point(i)));
+      value4 += l1.weight(i) * u.velocityNorm() * p2(u(l1.point(i)));
     }
 
-    REQUIRE(value == Catch::Approx(integrate(P, positions)));
+    REQUIRE(value == Catch::Approx(integrate(R, positions)));
+    REQUIRE(valuebis == Catch::Approx(integrate(Q, positions)));
     REQUIRE(value2 == Catch::Approx(integrate(P, positions2)));
     REQUIRE(value3 == Catch::Approx(integrate(P, positions3)));
     REQUIRE(value4 == Catch::Approx(integrate(P, positions4)));