diff --git a/maths/trapezoidal_rule.py b/maths/trapezoidal_rule.py
index 9a4ddc8af66b..3b872607943f 100644
--- a/maths/trapezoidal_rule.py
+++ b/maths/trapezoidal_rule.py
@@ -9,7 +9,7 @@
 """
 
 
-def method_1(boundary, steps):
+def method_1(func, boundary, steps):
     # "extended trapezoidal rule"
     # int(f) = dx/2 * (f1 + 2f2 + ... + fn)
     h = (boundary[1] - boundary[0]) / steps
@@ -20,21 +20,62 @@ def method_1(boundary, steps):
     y += (h / 2.0) * f(a)
     for i in x_i:
         # print(i)
-        y += h * f(i)
+        y += h * f(func, i)
     y += (h / 2.0) * f(b)
     return y
 
 
 def make_points(a, b, h):
-    x = a + h
-    while x < (b - h):
+    """
+    Yields equally spaced points between a and b
+    >>> x_i = make_points(0, 0.4, 0.1)
+    >>> next(x_i)
+    0.1
+    >>> next(x_i)
+    0.2
+    >>> next(x_i)
+    0.3
+    >>> next(x_i)
+    Traceback (most recent call last):
+        ...
+    StopIteration
+
+    >>> x_i = make_points(0, 1, -0.2)
+    >>> next(x_i)
+    Traceback (most recent call last):
+        ...
+    ValueError: h must be positive
+
+    >>> x_i = make_points(1,0,0.2)
+    >>> next(x_i)
+    Traceback (most recent call last):
+        ...
+    ValueError: a must be less than b
+    """
+    if h <= 0:
+        raise ValueError("h must be positive")
+    if a >= b:
+        raise ValueError("a must be less than b")
+
+    x = round(a + h, 10)
+    while x <= (b - h):
         yield x
-        x = x + h
+        x = round(x + h, 10)
 
 
-def f(x):  # enter your function here
-    y = (x - 0) * (x - 0)
-    return y
+def f(y, x):  # enter your function here
+    """
+    Returns the value of a lambda functiοn y at x
+    >>> f(lambda x: x**2,2)
+    4
+    >>> f(lambda x: x**2,-1)
+    1
+    >>> f(lambda x: (x+1)/2,5)
+    3.0
+    >>> f(lambda x: (x+1)/2,0)
+    0.5
+    """
+    return y(x)
 
 
 def main():
@@ -42,7 +83,7 @@ def main():
     b = 1.0  # Upper bound of integration
     steps = 10.0  # define number of steps or resolution
     boundary = [a, b]  # define boundary of integration
-    y = method_1(boundary, steps)
+    y = method_1(lambda x: x**2, boundary, steps)
     print(f"y = {y}")