diff --git a/DIRECTORY.md b/DIRECTORY.md
index 5117a9ba5..a5142cbec 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -163,6 +163,7 @@
   * [Gcd Iterative Euclidean](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/gcd_iterative_euclidean.cpp)
   * [Gcd Of N Numbers](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/gcd_of_n_numbers.cpp)
   * [Gcd Recursive Euclidean](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/gcd_recursive_euclidean.cpp)
+  * [Integral Approximation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/integral_approximation.cpp)
   * [Large Factorial](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/large_factorial.cpp)
   * [Large Number](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/large_number.h)
   * [Largest Power](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/largest_power.cpp)
diff --git a/math/integral_approximation.cpp b/math/integral_approximation.cpp
new file mode 100644
index 000000000..7e72ef404
--- /dev/null
+++ b/math/integral_approximation.cpp
@@ -0,0 +1,125 @@
+/**
+ * @file
+ * @brief Compute integral approximation of the function using [Riemann sum](https://en.wikipedia.org/wiki/Riemann_sum)
+ * @details In mathematics, a Riemann sum is a certain kind of approximation of an integral by a finite sum. It is named after nineteenth-century German mathematician Bernhard Riemann.
+ * One very common application is approximating the area of functions or lines on a graph and the length of curves and other approximations.
+ * The sum is calculated by partitioning the region into shapes (rectangles, trapezoids, parabolas, or cubics) that form a region similar to the region being measured, then calculating the area for each of these shapes, and finally adding all of these small areas together.
+ * This approach can be used to find a numerical approximation for a definite integral even if the fundamental theorem of calculus does not make it easy to find a closed-form solution.
+ * Because the region filled by the small shapes is usually not the same shape as the region being measured, the Riemann sum will differ from the area being measured.
+ * This error can be reduced by dividing up the region more finely, using smaller and smaller shapes. As the shapes get smaller and smaller, the sum approaches the Riemann integral.
+ * \author [Benjamin Walton](https://github.com/bwalton24)
+ * \author [Shiqi Sheng](https://github.com/shiqisheng00)
+ */
+#include <cassert>        /// for assert
+#include <cmath>         /// for mathematical functions
+#include <functional>   /// for passing in functions
+#include <iostream>    /// for IO operations
+
+/**
+ * @namespace math
+ * @brief Mathematical functions
+ */
+namespace math {
+/**
+ * @brief Computes integral approximation
+ * @param lb lower bound
+ * @param ub upper bound
+ * @param func function passed in
+ * @param delta
+ * @returns integral approximation of function from [lb, ub]
+ */
+double integral_approx(double lb, double ub,
+                       const std::function<double(double)>& func,
+                       double delta = .0001) {
+    double result = 0;
+    uint64_t numDeltas = static_cast<uint64_t>((ub - lb) / delta);
+    for (int i = 0; i < numDeltas; i++) {
+        double begin = lb + i * delta;
+        double end = lb + (i + 1) * delta;
+        result += delta * (func(begin) + func(end)) / 2;
+    }
+    return result;
+}
+
+/**
+ * @brief Wrapper to evaluate if the approximated
+ * value is within `.XX%` threshold of the exact value.
+ * @param approx aprroximate value
+ * @param exact expected value
+ * @param threshold values from [0, 1)
+ */
+void test_eval(double approx, double expected, double threshold) {
+    assert(approx >= expected * (1 - threshold));
+    assert(approx <= expected * (1 + threshold));
+}
+
+/**
+ * @brief Self-test implementations to
+ * test the `integral_approx` function.
+ *
+ * @returns `void`
+ */
+}  // namespace math
+
+static void test() {
+    double test_1 = math::integral_approx(
+        3.24, 7.56, [](const double x) { return log(x) + exp(x) + x; });
+    std::cout << "Test Case 1" << std::endl;
+    std::cout << "function: log(x) + e^x + x" << std::endl;
+    std::cout << "range: [3.24, 7.56]" << std::endl;
+    std::cout << "value: " << test_1 << std::endl;
+    math::test_eval(test_1, 1924.80384023549, .001);
+    std::cout << "Test 1 Passed!" << std::endl;
+    std::cout << "=====================" << std::endl;
+
+    double test_2 = math::integral_approx(0.023, 3.69, [](const double x) {
+        return x * x + cos(x) + exp(x) + log(x) * log(x);
+    });
+    std::cout << "Test Case 2" << std::endl;
+    std::cout << "function: x^2 + cos(x) + e^x + log^2(x)" << std::endl;
+    std::cout << "range: [.023, 3.69]" << std::endl;
+    std::cout << "value: " << test_2 << std::endl;
+    math::test_eval(test_2, 58.71291345202729, .001);
+    std::cout << "Test 2 Passed!" << std::endl;
+    std::cout << "=====================" << std::endl;
+
+    double test_3 = math::integral_approx(
+        10.78, 24.899, [](const double x) { return x * x * x - x * x + 378; });
+    std::cout << "Test Case 3" << std::endl;
+    std::cout << "function: x^3 - x^2 + 378" << std::endl;
+    std::cout << "range: [10.78, 24.899]" << std::endl;
+    std::cout << "value: " << test_3 << std::endl;
+    math::test_eval(test_3, 93320.65915078377, .001);
+    std::cout << "Test 3 Passed!" << std::endl;
+    std::cout << "=====================" << std::endl;
+
+    double test_4 = math::integral_approx(
+        .101, .505,
+        [](const double x) { return cos(x) * tan(x) * x * x + exp(x); },
+        .00001);
+    std::cout << "Test Case 4" << std::endl;
+    std::cout << "function: cos(x)*tan(x)*x^2 + e^x" << std::endl;
+    std::cout << "range: [.101, .505]" << std::endl;
+    std::cout << "value: " << test_4 << std::endl;
+    math::test_eval(test_4, 0.566485986311631, .001);
+    std::cout << "Test 4 Passed!" << std::endl;
+    std::cout << "=====================" << std::endl;
+
+    double test_5 = math::integral_approx(
+        -1, 1, [](const double x) { return exp(-1 / (x * x)); });
+    std::cout << "Test Case 5" << std::endl;
+    std::cout << "function: e^(-1/x^2)" << std::endl;
+    std::cout << "range: [-1, 1]" << std::endl;
+    std::cout << "value: " << test_5 << std::endl;
+    math::test_eval(test_5, 0.1781477117815607, .001);
+    std::cout << "Test 5 Passed!" << std::endl;
+}
+
+/**
+ * @brief Main function
+ * @returns 0 on exit
+ */
+int main() {
+    test();  // run self-test implementations
+    return 0;
+}