diff --git a/maths/primelib.py b/maths/primelib.py
index cf01750cf912..bb5c21078f5b 100644
--- a/maths/primelib.py
+++ b/maths/primelib.py
@@ -38,7 +38,6 @@
 """
 
 from math import sqrt
-
 from maths.greatest_common_divisor import gcd_by_iterative
 
 
@@ -46,6 +45,15 @@ def is_prime(number: int) -> bool:
     """
     input: positive integer 'number'
     returns true if 'number' is prime otherwise false.
+
+    >>> is_prime(97)
+    True
+    >>> is_prime(9991)
+    False
+    >>> is_prime(-3)
+    Traceback (most recent call last):
+        ...
+    AssertionError: 'number' must been an int and positive
     """
 
     # precondition
@@ -589,3 +597,9 @@ def fib(n):
         fib1 = tmp
 
     return ans
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/maths/stewart.py b/maths/stewart.py
new file mode 100644
index 000000000000..edae5f61e959
--- /dev/null
+++ b/maths/stewart.py
@@ -0,0 +1,47 @@
+"""
+In geometry, Stewart's theorem yields a relation between the
+lengths of the sides and the length of a cevian in a triangle.
+Its name is in honour of the Scottish mathematician Matthew
+Stewart, who published the theorem in 1746.[1]
+
+Source: https://en.wikipedia.org/wiki/Stewart%27s_theorem
+"""
+
+
+def stewart(lena: float, lenb: float, lenc: float, lenn: float, lenm: float) -> float:
+    """
+    Given the side lengths of the triangle (a,b,c), where the cevian intersects
+    the side with side length a and splits it into segments with lengths n and m,
+    this formula finds the length of the cevian.
+
+    >>> stewart(1,1,1,0.5,0.5)
+    0.8660254037844386
+    >>> stewart(1,2,3,4,5)
+    Traceback (most recent call last):
+        ...
+    ValueError: This triangle violates the triangle inequality
+    >>> stewart(1,1,1,1,1)
+    Traceback (most recent call last):
+        ...
+    ValueError: n+m must equal a
+    >>> stewart(3,2,4,1.7,1.3)
+    2.9308701779505686
+    """
+    a = lena
+    b = lenb
+    c = lenc
+    n = lenn
+    m = lenm
+    if a + b <= c or b + c <= a or a + c <= b:
+        raise ValueError("This triangle violates the triangle inequality")
+    if n + m != a:
+        raise ValueError("n+m must equal a")
+    if a <= 0 or b <= 0 or c <= 0 or n < 0 or m < 0:
+        raise ValueError("The side lengths of a triangle have to be positive")
+    return ((b**2 * m + c**2 * n - m * a * n) / a) ** 0.5
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()