diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs
index 84398dbe..01400e98 100644
--- a/src/Numerics/Integrate.cs
+++ b/src/Numerics/Integrate.cs
@@ -62,5 +62,33 @@ namespace MathNet.Numerics
         {
             return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, 1e-8);
         }
+
+        /// <summary>
+        /// Approximates a definite integral using an Nth order Gauss-Legendre rule.
+        /// </summary>
+        /// <param name="f">The analytic smooth function to integrate.</param>
+        /// <param name="invervalBegin">Where the interval starts, exclusive and finite.</param>
+        /// <param name="invervalEnd">Where the interval ends, exclusive and finite.</param>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2,. . ., 20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calulcated on the fly.</param>
+        /// <returns>Approximation of the finite integral in the given interval.</returns>
+        public static double GaussLegendre(Func<double, double> f, double invervalBegin, double invervalEnd, int order)
+        {
+            return GaussLegendreRule.Integrate(f, invervalBegin, invervalEnd, order);
+        }
+
+        /// <summary>
+        /// Approximates a 2-dimensional definite integral using an Nth order Gauss-Legendre rule over the rectangle [a,b] x [c,d].
+        /// </summary>
+        /// <param name="f">The 2-dimensional analytic smooth function to integrate.</param>
+        /// <param name="invervalBeginA">Where the interval starts for the first (inside) integral, exclusive and finite.</param>
+        /// <param name="invervalEndA">Where the interval ends for the first (inside) integral, exclusive and finite.</param>
+        /// <param name="invervalBeginB">Where the interval starts for the second (outside) integral, exclusive and finite.</param>
+        /// /// <param name="invervalEndB">Where the interval ends for the second (outside) integral, exclusive and finite.</param>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule. Precomputed Gauss-Legendre abscissas/weights for orders 2,. . ., 20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calulcated on the fly.</param>
+        /// <returns>Approximation of the finite integral in the given interval.</returns>
+        public static double GaussLegendre(Func<double, double, double> f, double invervalBeginA, double invervalEndA, double invervalBeginB, double invervalEndB, int order)
+        {
+            return GaussLegendreRule.Integrate(f, invervalBeginA, invervalEndA, invervalBeginB, invervalEndB, order);
+        }
     }
 }
diff --git a/src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs b/src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs
index 78967b4b..3dbd961f 100644
--- a/src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs
+++ b/src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs
@@ -44,7 +44,7 @@ namespace MathNet.Numerics.Integration.GaussRule
         /// <summary>
         /// Getter for the GaussPoint.
         /// </summary>
-        /// <param name="order">Defines an Nth order Gauss-Legendre rule. Precomputed Gauss-Legendre abscissas/weights for orders 2,. . ., 20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calulcated on the fly. Furthermore, caching in a private static field is used to speed up the algorithm.</param>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. Precomputed Gauss-Legendre abscissas/weights for orders 2,. . ., 20, 32, 64, 96, 100, 128, 256, 512, 1024 are used, otherwise they're calulcated on the fly.</param>
         /// <returns>Object containing the non-negative abscissas/weights, order, and intervalBegin/intervalEnd. The non-negative abscissas/weights are generated over the interval [-1,1] for the given order.</returns>
         public static GaussPoint GetGaussPoint(int order)
         {
diff --git a/src/UnitTests/IntegrationTests/IntegrationTest.cs b/src/UnitTests/IntegrationTests/IntegrationTest.cs
index 26c6851c..2a6f3fe0 100644
--- a/src/UnitTests/IntegrationTests/IntegrationTest.cs
+++ b/src/UnitTests/IntegrationTests/IntegrationTest.cs
@@ -220,6 +220,21 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
             Assert.Less(relativeError, 5e-16);
         }
 
+        /// <summary>
+        /// Gauss-Legendre rule supports integration.
+        /// </summary>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</param>
+        [TestCase(19)]
+        [TestCase(20)]
+        [TestCase(21)]
+        [TestCase(22)]
+        public void TestIntegrateGaussLegendre(int order)
+        {
+            double appoximateArea = Integrate.GaussLegendre(TargetFunctionA, StartA, StopA, order);
+            double relativeError = Math.Abs(TargetAreaA - appoximateArea) / TargetAreaA;
+            Assert.Less(relativeError, 5e-16);
+        }
+
         /// <summary>
         /// Gauss-Legendre rule supports 2-dimensional integration over the rectangle.
         /// </summary>
@@ -235,6 +250,21 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
             Assert.Less(relativeError, 1e-15);
         }
 
+        /// <summary>
+        /// Gauss-Legendre rule supports 2-dimensional integration over the rectangle.
+        /// </summary>
+        /// <param name="order">Defines an Nth order Gauss-Legendre rule. The order also defines the number of abscissas and weights for the rule.</param>
+        [TestCase(19)]
+        [TestCase(20)]
+        [TestCase(21)]
+        [TestCase(22)]
+        public void TestIntegrateGaussLegendre2D(int order)
+        {
+            double appoximateArea = Integrate.GaussLegendre(TargetFunctionB, StartA, StopA, StartB, StopB, order);
+            double relativeError = Math.Abs(TargetAreaB - appoximateArea) / TargetAreaB;
+            Assert.Less(relativeError, 1e-15);
+        }
+
         /// <summary>
         /// Gauss-Legendre rule supports obtaining the abscissas/weights. In this case, they're used for integration. 
         /// </summary>