diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 3ef5e180710d90e0b62e0cf8f8d222cb3a84ca14..c6cf4f8165e2f7a4eb812b85d0ee4097b197169c 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -464,4 +464,81 @@ inverse(const TinyMatrix<3, T>& A)
   return A_cofactors_T *= 1. / determinent;
 }
 
+PUGS_INLINE constexpr void
+_swap(double& a, double& b)
+{
+  double temp = a;
+  a           = b;
+  b           = temp;
+}
+
+template <size_t N, typename T>
+PUGS_INLINE constexpr TinyMatrix<N, T>
+inverse(const TinyMatrix<N, T>& A)
+{
+  static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
+  static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
+  TinyVector<N, size_t> indexc(zero);
+  TinyVector<N, size_t> indexr(zero);
+  TinyVector<N, size_t> indpiv(zero);
+  TinyMatrix<N> inv_A(A);
+  double big    = 0;
+  double pivinv = 0;
+  double dum    = 0;
+  size_t irow   = 0;
+  size_t icol   = 0;
+
+  for (size_t i = 0; i < N; ++i) {
+    big = 0;
+    for (size_t j = 0; j < N; ++j) {
+      if (indpiv[j] != 1) {
+        for (size_t k = 0; k < N; k++) {
+          if (indpiv[k] == 0) {
+            if (std::abs(inv_A(j, k)) >= big) {
+              big  = std::abs(inv_A(j, k));
+              irow = j;
+              icol = k;
+            }
+          } else {
+            Assert(indpiv[k] <= 1, "Gauss pivoting: singular matrix");
+          }
+        }
+      }
+    }
+    ++indpiv[icol];
+    // We found the pivot A(irow,icol), now we swap columns to put the pivot on the diagonal
+    if (irow != icol) {
+      for (size_t l = 0; l < N; ++l) {
+        _swap(inv_A(irow, l), inv_A(icol, l));
+      }
+    }
+    // we save the swap to invert it at the end
+    indexr[i] = irow;
+    indexc[i] = icol;
+    Assert(inv_A(icol, icol) != 0, "Gauss pivoting: singular matrix");
+    pivinv            = 1 / inv_A(icol, icol);
+    inv_A(icol, icol) = 1;
+    for (size_t l = 0; l < N; ++l) {
+      inv_A(icol, l) *= pivinv;
+    }
+    for (size_t ll = 0; ll < N; ++ll) {
+      if (ll != icol) {
+        dum             = inv_A(ll, icol);
+        inv_A(ll, icol) = 0;
+        for (size_t l = 0; l < N; ++l) {
+          inv_A(ll, l) -= inv_A(icol, l) * dum;
+        }
+      }
+    }
+  }
+  for (size_t l = N; l > 0; --l) {
+    if (indexr[l - 1] != indexc[l - 1]) {
+      for (size_t k = 0; k < N; ++k) {
+        _swap(inv_A(k, indexr[l - 1]), inv_A(k, indexc[l - 1]));
+      }
+    }
+  }
+  return inv_A;
+}
+
 #endif   // TINYMATRIX_HPP
diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp
index 653bf69f8bea6af7b242e954b48698160d143572..6a2bce5ec081ce80c40e7dfc1c4598f7abd70b2a 100644
--- a/src/analysis/Polynomial.hpp
+++ b/src/analysis/Polynomial.hpp
@@ -328,29 +328,6 @@ class Polynomial
     return Polynomial<N + 1>{coefs};
   }
 
-  PUGS_INLINE
-  constexpr friend double
-  integrate(const Polynomial<N>& P, const double& xinf, const double& xsup)
-  {
-    Polynomial<N + 1> Q = primitive(P);
-    return (Q(xsup) - Q(xinf));
-  }
-
-  PUGS_INLINE
-  constexpr friend auto
-  derivative(const Polynomial<N>& P)
-  {
-    if constexpr (N == 0) {
-      return Polynomial<0>(0);
-    } else {
-      TinyVector<N> coefs;
-      for (size_t i = 0; i < N; ++i) {
-        coefs[i] = double(i + 1) * P.coefficients()[i + 1];
-      }
-      return Polynomial<N - 1>(coefs);
-    }
-  }
-
   PUGS_INLINE
   constexpr friend std::ostream&
   operator<<(std::ostream& os, const Polynomial<N>& P)
@@ -410,31 +387,13 @@ class Polynomial
     return os;
   }
 
-  PUGS_INLINE
-  constexpr friend void
-  lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k, Polynomial<N>& lk)
-  {
-    for (size_t i = 0; i < N; ++i) {
-      Assert(zeros[i] < zeros[i + 1], "Interpolation values must be strictly increasing in Lagrange polynomials");
-    }
-    lk.coefficients()    = zero;
-    lk.coefficients()[0] = 1;
-    for (size_t i = 0; i < N + 1; ++i) {
-      if (k == i)
-        continue;
-      double factor = 1. / (zeros[k] - zeros[i]);
-      Polynomial<1> P({-zeros[i] * factor, factor});
-      lk *= P;
-    }
-  }
-
   PUGS_INLINE
   constexpr friend void
   lagrangeBasis(const TinyVector<N + 1> zeros, TinyVector<N + 1, Polynomial<N>>& basis)
   {
     Polynomial<N> lj;
     for (size_t j = 0; j < N + 1; ++j) {
-      lagrangePolynomial(zeros, j, basis[j]);
+      basis[j] = lagrangePolynomial(zeros, j);
     }
   }
 
@@ -460,6 +419,9 @@ class Polynomial
   ~Polynomial()                   = default;
 };
 
+template <size_t N>
+PUGS_INLINE constexpr Polynomial<N> lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k);
+
 template <size_t N>
 PUGS_INLINE constexpr TinyVector<N, Polynomial<N - 1>>
 lagrangeBasis(const TinyVector<N>& zeros)
@@ -468,9 +430,65 @@ lagrangeBasis(const TinyVector<N>& zeros)
   TinyVector<N, Polynomial<N - 1>> basis;
   Polynomial<N - 1> lj;
   for (size_t j = 0; j < N; ++j) {
-    lagrangePolynomial(zeros, j, basis[j]);
+    basis[j] = lagrangePolynomial<N - 1>(zeros, j);
   }
   return basis;
 }
 
+template <size_t N>
+PUGS_INLINE constexpr double
+integrate(const Polynomial<N>& P, const double& xinf, const double& xsup)
+{
+  Polynomial<N + 1> Q = primitive(P);
+  return (Q(xsup) - Q(xinf));
+}
+
+template <size_t N>
+PUGS_INLINE constexpr double
+symmetricIntegrate(const Polynomial<N>& P, const double& delta)
+{
+  Assert(delta > 0, "A positive delta is needed for symmetricIntegrate");
+  double integral = 0.;
+  for (size_t i = 0; i <= N; ++i) {
+    if (i % 2 == 0)
+      integral += 2. * P.coefficients()[i] * std::pow(delta, i + 1) / (i + 1);
+  }
+  return integral;
+}
+
+template <size_t N>
+PUGS_INLINE constexpr auto
+derivative(const Polynomial<N>& P)
+{
+  if constexpr (N == 0) {
+    return Polynomial<0>(0);
+  } else {
+    TinyVector<N> coefs;
+    for (size_t i = 0; i < N; ++i) {
+      coefs[i] = double(i + 1) * P.coefficients()[i + 1];
+    }
+    return Polynomial<N - 1>(coefs);
+  }
+}
+
+template <size_t N>
+PUGS_INLINE constexpr Polynomial<N>
+lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k)
+{
+  for (size_t i = 0; i < N; ++i) {
+    Assert(zeros[i] < zeros[i + 1], "Interpolation values must be strictly increasing in Lagrange polynomials");
+  }
+  Polynomial<N> lk;
+  lk.coefficients()    = zero;
+  lk.coefficients()[0] = 1;
+  for (size_t i = 0; i < N + 1; ++i) {
+    if (k == i)
+      continue;
+    double factor = 1. / (zeros[k] - zeros[i]);
+    Polynomial<1> P({-zeros[i] * factor, factor});
+    lk *= P;
+  }
+  return lk;
+}
+
 #endif   // POLYNOMIAL_HPP
diff --git a/src/analysis/PolynomialBasis.hpp b/src/analysis/PolynomialBasis.hpp
index c8e782b799a8b6baecbaf80af594d9d72ca6ee5a..4359f408d93db2bc814a9225f69b73982d04f796 100644
--- a/src/analysis/PolynomialBasis.hpp
+++ b/src/analysis/PolynomialBasis.hpp
@@ -1,5 +1,5 @@
-#ifndef POLYNOMIALBASIS_HPP
-#define POLYNOMIALBASIS_HPP
+#ifndef POLYNOMIAL_BASIS_HPP
+#define POLYNOMIAL_BASIS_HPP
 
 #include <algebra/TinyVector.hpp>
 #include <analysis/Polynomial.hpp>
@@ -7,7 +7,9 @@
 
 enum class BasisType
 {
-  canonical
+  lagrange,
+  canonical,
+  taylor
 };
 
 template <size_t N>
@@ -29,6 +31,46 @@ class PolynomialBasis
     return *this;
   }
 
+  PUGS_INLINE
+  constexpr PolynomialBasis<N>&
+  _buildTaylorBasis(const double& shift)
+  {
+    TinyVector<N + 1> coefficients(zero);
+    elements()[0] = Polynomial<N>(coefficients);
+    elements()[0] += Polynomial<0>(TinyVector<1>{1});
+    Polynomial<1> unit(TinyVector<2>{-shift, 1});
+    for (size_t i = 1; i <= N; i++) {
+      elements()[i] = elements()[i - 1] * unit;
+    }
+    return *this;
+  }
+
+  PUGS_INLINE
+  constexpr PolynomialBasis<N>&
+  _buildLagrangeBasis(const TinyVector<N + 1>& zeros)
+  {
+    for (size_t i = 0; i < N; ++i) {
+      Assert(zeros[i] < zeros[i + 1], "Interpolation values must be strictly increasing in Lagrange polynomials");
+    }
+    if constexpr (N == 0) {
+      elements()[0] = Polynomial<0>(TinyVector<1>{1});
+      return *this;
+    } else {
+      for (size_t i = 0; i <= N; ++i) {
+        TinyVector<N + 1> coefficients(zero);
+        elements()[i]                   = Polynomial<N>(coefficients);
+        elements()[i].coefficients()[0] = 1;
+        for (size_t j = 0; j < N + 1; ++j) {
+          if (i == j)
+            continue;
+          double adim = 1. / (zeros[i] - zeros[j]);
+          elements()[i] *= Polynomial<1>{{-zeros[j] * adim, adim}};
+        }
+      }
+      return *this;
+    }
+  }
+
  public:
   PUGS_INLINE
   constexpr size_t
@@ -56,8 +98,12 @@ class PolynomialBasis
   displayType()
   {
     switch (m_basis_type) {
+    case BasisType::lagrange:
+      return "lagrange";
     case BasisType::canonical:
       return "canonical";
+    case BasisType::taylor:
+      return "taylor";
     default:
       return "unknown basis type";
     }
@@ -79,14 +125,22 @@ class PolynomialBasis
 
   PUGS_INLINE
   constexpr PolynomialBasis<N>&
-  build(BasisType basis_type)
+  build(BasisType basis_type, const double& shift = 0, const TinyVector<N + 1>& zeros = TinyVector<N + 1>(zero))
   {
     type() = basis_type;
     switch (basis_type) {
+    case BasisType::lagrange: {
+      return _buildLagrangeBasis(zeros);
+      break;
+    }
     case BasisType::canonical: {
       return _buildCanonicalBasis();
       break;
     }
+    case BasisType::taylor: {
+      return _buildTaylorBasis(shift);
+      break;
+    }
     default:
       throw UnexpectedError("unknown basis type");
     }
@@ -102,4 +156,4 @@ class PolynomialBasis
   constexpr PolynomialBasis() noexcept = default;
   ~PolynomialBasis()                   = default;
 };
-#endif   // POLYNOMIAL_HPP
+#endif   // POLYNOMIAL_BASIS_HPP
diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
index a8599f575058185b5ceeed0814b242ae6ec70b0a..52dfba5dd4576f4739d44dc6004601e144a88b0d 100644
--- a/src/language/algorithms/CMakeLists.txt
+++ b/src/language/algorithms/CMakeLists.txt
@@ -2,6 +2,7 @@
 
 add_library(PugsLanguageAlgorithms
   AcousticSolverAlgorithm.cpp
+  ODEDiscontinuousGalerkin1D.cpp
 )
 
 
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 542690a4ffcfb73ba174f0014b358548208d830f..a86a752c2dbb3d7be4304e638a3f7c269311bcb1 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,6 +1,7 @@
 #include <language/modules/SchemeModule.hpp>
 
 #include <language/algorithms/AcousticSolverAlgorithm.hpp>
+#include <language/algorithms/ODEDiscontinuousGalerkin1D.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Mesh.hpp>
@@ -173,5 +174,60 @@ SchemeModule::SchemeModule()
                                 }
                               }
 
+                              ));
+
+  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const BasisType>>);
+
+  this->_addBuiltinFunction("canonicalBasis",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const BasisType>()>>(
+
+                              []() -> std::shared_ptr<BasisType> {
+                                return std::make_shared<BasisType>(BasisType::canonical);
+                              }
+
+                              ));
+  this->_addBuiltinFunction("taylorBasis",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const BasisType>()>>(
+
+                              []() -> std::shared_ptr<BasisType> {
+                                return std::make_shared<BasisType>(BasisType::taylor);
+                              }
+
+                              ));
+  this->_addBuiltinFunction("lagrangeBasis",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const BasisType>()>>(
+
+                              []() -> std::shared_ptr<BasisType> {
+                                return std::make_shared<BasisType>(BasisType::lagrange);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("odediscontinuousgalerkin1D",
+                            std::make_shared<
+                              BuiltinFunctionEmbedder<void(std::shared_ptr<const BasisType> basis_type, const size_t&,
+                                                           std::shared_ptr<const IMesh>, const FunctionSymbolId&)>>(
+                              [](std::shared_ptr<const BasisType> basis_type, const size_t& Degree,
+                                 std::shared_ptr<const IMesh> p_mesh, const FunctionSymbolId& uex_id) -> void {
+                                Assert(p_mesh->dimension() == 1, "invalid mesh dimension");
+                                switch (Degree) {
+                                case 0: {
+                                  ODEDiscontinuousGalerkin1D<1, 0>(*basis_type, p_mesh, uex_id);
+                                  break;
+                                }
+                                case 1: {
+                                  ODEDiscontinuousGalerkin1D<1, 1>(*basis_type, p_mesh, uex_id);
+                                  break;
+                                }
+                                case 2: {
+                                  ODEDiscontinuousGalerkin1D<1, 2>(*basis_type, p_mesh, uex_id);
+                                  break;
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid polynomial degree");
+                                }
+                                }
+                              }
+
                               ));
 }
diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp
index 3c70be9a2db29d634b71ea96fb60fce38183c85d..b140295b486a33128295fb9829d2c6d540277948 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -1,6 +1,7 @@
 #ifndef SCHEME_MODULE_HPP
 #define SCHEME_MODULE_HPP
 
+#include <analysis/PolynomialBasis.hpp>
 #include <language/modules/BuiltinModule.hpp>
 #include <language/utils/ASTNodeDataTypeTraits.hpp>
 #include <utils/PugsMacros.hpp>
@@ -15,6 +16,10 @@ template <>
 inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>> =
   {ASTNodeDataType::type_id_t, "boundary_condition"};
 
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const BasisType>> = {ASTNodeDataType::type_id_t,
+                                                                                    "basis_type"};
+
 class SchemeModule : public BuiltinModule
 {
  public:
diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp
index fd1f07a2c24cb8113de96feb5049f94682e3f266..5e9eb27bb41609b9da10fb81b79a90c19fd49830 100644
--- a/tests/test_Polynomial.cpp
+++ b/tests/test_Polynomial.cpp
@@ -126,6 +126,8 @@ TEST_CASE("Polynomial", "[analysis]")
     double xsup   = 1;
     double result = integrate(P, xinf, xsup);
     REQUIRE(result == 6);
+    result = symmetricIntegrate(P, 2);
+    REQUIRE(result == 24);
   }
   SECTION("derivative")
   {
@@ -182,17 +184,13 @@ TEST_CASE("Polynomial", "[analysis]")
   {
     Polynomial<1> S({0.5, -0.5});
     Polynomial<1> Q;
-    lagrangePolynomial({-1, 1}, 0, Q);
+    Q = lagrangePolynomial<1>({-1, 1}, 0);
     REQUIRE(S == Q);
     Polynomial<2> P({0, -0.5, 0.5});
     Polynomial<2> R;
-    lagrangePolynomial({-1, 0, 1}, 0, R);
+    R = lagrangePolynomial<2>({-1, 0, 1}, 0);
     REQUIRE(R == P);
     const TinyVector<3, Polynomial<2>> basis = lagrangeBasis(TinyVector<3>{-1, 0, 1});
-    //    lagrangeBasis({-1, 0, 1}, basis);
-    // basis[0] = Polynomial<2>({0, -0.5, 0.5});
-    // basis[1] = Polynomial<2>({1, 0, -1});
-    // basis[2] = Polynomial<2>({0, 0.5, 0.5});
     REQUIRE(lagrangeToCanonical({1, 0, 1}, basis) == Polynomial<2>({0, 0, 1}));
   }
 }
diff --git a/tests/test_PolynomialBasis.cpp b/tests/test_PolynomialBasis.cpp
index 933797bb328386809133fe5c7d14c5cdb318a8c9..ae3745553d9e5dc4a528e32a86818bd7965301e0 100644
--- a/tests/test_PolynomialBasis.cpp
+++ b/tests/test_PolynomialBasis.cpp
@@ -38,5 +38,10 @@ TEST_CASE("PolynomialBasis", "[analysis]")
     B.build(BasisType::canonical);
     REQUIRE(B.elements()[1] == Polynomial<2>{{0, 1, 0}});
     REQUIRE(B.elements()[2] == Polynomial<2>{{0, 0, 1}});
+    PolynomialBasis<2> C;
+    C.build(BasisType::taylor);
+    REQUIRE(B.elements()[1] == C.elements()[1]);
+    C.build(BasisType::taylor, 1);
+    REQUIRE(C.elements()[2] == Polynomial<2>{{1, -2, 1}});
   }
 }
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index e0d7ef39da53ae03358b307e448136d69216d1cf..e3245e79612532f0919d142fe0a349ce69828b4d 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -168,7 +168,7 @@ TEST_CASE("TinyMatrix", "[algebra]")
 
   SECTION("checking for inverse calculations")
   {
-    {
+    {   // // //
       const TinyMatrix<1, double> A1(2);
       REQUIRE(inverse(A1)(0, 0) == Approx(0.5).epsilon(1E-14));
     }
@@ -197,6 +197,27 @@ TEST_CASE("TinyMatrix", "[algebra]")
       REQUIRE(I(2, 1) == Approx(0).margin(1E-14));
       REQUIRE(I(2, 2) == Approx(1).epsilon(1E-14));
     }
+    {
+      const TinyMatrix<4, double> A4(2, 3, 1, 4, -1, 5, -2, 3, 4, 1, 2, 3, 4, -1, -1, 0);
+      const TinyMatrix<4, double> I = inverse(A4) * A4;
+
+      REQUIRE(I(0, 0) == Approx(1).epsilon(1E-14));
+      REQUIRE(I(0, 1) == Approx(0).margin(1E-14));
+      REQUIRE(I(0, 2) == Approx(0).margin(1E-14));
+      REQUIRE(I(0, 3) == Approx(0).margin(1E-14));
+      REQUIRE(I(1, 0) == Approx(0).margin(1E-14));
+      REQUIRE(I(1, 1) == Approx(1).epsilon(1E-14));
+      REQUIRE(I(1, 2) == Approx(0).margin(1E-14));
+      REQUIRE(I(1, 3) == Approx(0).margin(1E-14));
+      REQUIRE(I(2, 0) == Approx(0).margin(1E-14));
+      REQUIRE(I(2, 1) == Approx(0).margin(1E-14));
+      REQUIRE(I(2, 2) == Approx(1).epsilon(1E-14));
+      REQUIRE(I(2, 3) == Approx(0).margin(1E-14));
+      REQUIRE(I(3, 0) == Approx(0).margin(1E-14));
+      REQUIRE(I(3, 1) == Approx(0).margin(1E-14));
+      REQUIRE(I(3, 2) == Approx(0).margin(1E-14));
+      REQUIRE(I(3, 3) == Approx(1).epsilon(1E-14));
+    }
   }
 
   SECTION("checking for matrices output")