diff --git a/DIRECTORY.md b/DIRECTORY.md
index f0a34a553946..02b149525484 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -561,6 +561,7 @@
 ## Linear Algebra
   * [Gaussian Elimination](linear_algebra/gaussian_elimination.py)
   * [Jacobi Iteration Method](linear_algebra/jacobi_iteration_method.py)
+  * [Lanczos Algorithm](linear_algebra/lanczos_algorithm.py)
   * [Lu Decomposition](linear_algebra/lu_decomposition.py)
   * Src
     * [Conjugate Gradient](linear_algebra/src/conjugate_gradient.py)
@@ -719,6 +720,7 @@
   * [Pollard Rho](maths/pollard_rho.py)
   * [Polynomial Evaluation](maths/polynomial_evaluation.py)
   * Polynomials
+    * [Legendre](maths/polynomials/legendre.py)
     * [Single Indeterminate Operations](maths/polynomials/single_indeterminate_operations.py)
   * [Power Using Recursion](maths/power_using_recursion.py)
   * [Prime Check](maths/prime_check.py)
diff --git a/linear_algebra/lanczos_algorithm.py b/linear_algebra/lanczos_algorithm.py
new file mode 100644
index 000000000000..f59986da756e
--- /dev/null
+++ b/linear_algebra/lanczos_algorithm.py
@@ -0,0 +1,37 @@
+import numpy as np
+
+
+def lanczos(a: np.ndarray) -> tuple[list[float], list[float]]:
+    """
+    Implements the Lanczos algorithm for a symmetric matrix.
+
+    Parameters:
+    -----------
+    matrix : numpy.ndarray
+        Symmetric matrix of size (n, n).
+
+    Returns:
+    --------
+    alpha : [float]
+        List of diagonal elements of the resulting tridiagonal matrix.
+    beta : [float]
+        List of off-diagonal elements of the resulting tridiagonal matrix.
+    """
+    n = a.shape[0]
+    v = np.zeros((n, n))
+    rng = np.random.default_rng()
+    v[:, 0] = rng.standard_normal(n)
+    v[:, 0] /= np.linalg.norm(v[:, 0])
+    alpha: list[float] = []
+    beta: list[float] = []
+    for j in range(n):
+        w = np.dot(a, v[:, j])
+        alpha.append(np.dot(w, v[:, j]))
+        if j == n - 1:
+            break
+        w -= alpha[j] * v[:, j]
+        if j > 0:
+            w -= beta[j - 1] * v[:, j - 1]
+        beta.append(np.linalg.norm(w))
+        v[:, j + 1] = w / beta[j]
+    return alpha, beta
diff --git a/maths/polynomials/legendre.py b/maths/polynomials/legendre.py
new file mode 100644
index 000000000000..8bd3fdb56832
--- /dev/null
+++ b/maths/polynomials/legendre.py
@@ -0,0 +1,60 @@
+from numpy.polynomial import Polynomial
+from math import factorial
+import pytest
+
+
+def legendre(n: int) -> list[float]:
+    """
+    Compute the coefficients of the nth Legendre polynomial.
+
+    The Legendre polynomials are solutions to Legendre's differential equation
+    and are widely used in physics and engineering.
+
+    Parameters:
+        n (int): The order of the Legendre polynomial.
+
+    Returns:
+        list[float]: Coefficients of the polynomial in ascending order of powers.
+    """
+    p = (1 / (factorial(n) * (2 ** n))) * (Polynomial([-1, 0, 1]) ** n)
+    return p.deriv(n).coef.tolist()
+
+
+def jsp():
+    print(legendre(1))
+    print(legendre(2))
+    print(legendre(3))
+    print(legendre(4))
+
+
+jsp()
+
+
+def test_legendre_0():
+    """Test the 0th Legendre polynomial."""
+    assert legendre(0) == [1.0], "The 0th Legendre polynomial should be [1.0]"
+
+
+def test_legendre_1():
+    """Test the 1st Legendre polynomial."""
+    assert legendre(1) == [0.0, 1.0], "The 1st Legendre polynomial should be [0.0, 1.0]"
+
+
+def test_legendre_2():
+    """Test the 2nd Legendre polynomial."""
+    assert legendre(2) == [-0.5, 0.0, 1.5], "The 2nd Legendre polynomial should be [-0.5, 0.0, 1.5]"
+
+
+def test_legendre_3():
+    """Test the 3rd Legendre polynomial."""
+    assert legendre(3) == [0.0, -1.5, 0.0, 2.5], "The 3rd Legendre polynomial should be [0.0, -1.5, 0.0, 2.5]"
+
+
+def test_legendre_4():
+    """Test the 4th Legendre polynomial."""
+    assert legendre(4) == pytest.approx([0.375, 0.0, -3.75, 0.0, 4.375])
+    "The 4th Legendre polynomial should be [0.375, 0.0, -3.75, 0.0, 4.375]"
+
+
+if __name__ == '__main__':
+    pytest.main()