diff --git a/integral/simpson_integral.py b/integral/simpson_integral.py
new file mode 100644
index 0000000..2eafac6
--- /dev/null
+++ b/integral/simpson_integral.py
@@ -0,0 +1,39 @@
+import numpy as np
+
+def simpson_rule_1_3(a, b, f, n=100):
+    """
+    Function to calculation simpson 1/3
+    >>> f = lambda x: x**3
+    >>> simpson_rule_1_3(0, 2, f, n=100)
+    4.000000000000002
+    """
+    res = 0
+    h = (b - a)/ n
+    for i in range(0, n, 2):
+        x0 = a + i*h
+        x1 = x0 + h
+        x2 = x0 + 2 * h
+        res += (h/3) * (f(x0) + 4*f(x1) + f(x2))
+    return res
+
+def simpson_rule_3_8(a, b, f, n=100):
+    """
+    Function to calculation simpson 3/8
+    >>> f = lambda x: x**2
+    >>> simpson_rule_3_8(0, 1, f, n=99)
+    0.3333333333333334
+    """
+    h = (b - a) / n
+    res = 0
+    for i in range(0, n, 3):
+        x0 = a + i*h
+        x1 = x0 + h
+        x2 = x0 + 2*h
+        x3 = x0 + 3*h
+        res += (3*h/8) * (f(x0) + 3*f(x1) + 3*f(x2) + f(x3))
+    return res
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()