diff --git a/linear_algebra/gaussian_elimination.py b/linear_algebra/gaussian_elimination.py
index 724773c0db98..6f4075b710fd 100644
--- a/linear_algebra/gaussian_elimination.py
+++ b/linear_algebra/gaussian_elimination.py
@@ -1,6 +1,6 @@
 """
-Gaussian elimination method for solving a system of linear equations.
-Gaussian elimination - https://en.wikipedia.org/wiki/Gaussian_elimination
+| Gaussian elimination method for solving a system of linear equations.
+| Gaussian elimination - https://en.wikipedia.org/wiki/Gaussian_elimination
 """
 
 import numpy as np
@@ -13,12 +13,17 @@ def retroactive_resolution(
 ) -> NDArray[float64]:
     """
     This function performs a retroactive linear system resolution
-        for triangular matrix
+    for triangular matrix
 
     Examples:
-        2x1 + 2x2 - 1x3 = 5         2x1 + 2x2 = -1
-        0x1 - 2x2 - 1x3 = -7        0x1 - 2x2 = -1
-        0x1 + 0x2 + 5x3 = 15
+        1.
+            * 2x1 + 2x2 - 1x3 = 5
+            * 0x1 - 2x2 - 1x3 = -7
+            * 0x1 + 0x2 + 5x3 = 15
+        2.
+            * 2x1 + 2x2 = -1
+            * 0x1 - 2x2 = -1
+
     >>> gaussian_elimination([[2, 2, -1], [0, -2, -1], [0, 0, 5]], [[5], [-7], [15]])
     array([[2.],
            [2.],
@@ -45,9 +50,14 @@ def gaussian_elimination(
     This function performs Gaussian elimination method
 
     Examples:
-        1x1 - 4x2 - 2x3 = -2        1x1 + 2x2 = 5
-        5x1 + 2x2 - 2x3 = -3        5x1 + 2x2 = 5
-        1x1 - 1x2 + 0x3 = 4
+        1.
+            * 1x1 - 4x2 - 2x3 = -2
+            * 5x1 + 2x2 - 2x3 = -3
+            * 1x1 - 1x2 + 0x3 = 4
+        2.
+            * 1x1 + 2x2 = 5
+            * 5x1 + 2x2 = 5
+
     >>> gaussian_elimination([[1, -4, -2], [5, 2, -2], [1, -1, 0]], [[-2], [-3], [4]])
     array([[ 2.3 ],
            [-1.7 ],
diff --git a/linear_algebra/lu_decomposition.py b/linear_algebra/lu_decomposition.py
index 3620674835cd..3d89b53a48fb 100644
--- a/linear_algebra/lu_decomposition.py
+++ b/linear_algebra/lu_decomposition.py
@@ -2,13 +2,14 @@
 Lower-upper (LU) decomposition factors a matrix as a product of a lower
 triangular matrix and an upper triangular matrix. A square matrix has an LU
 decomposition under the following conditions:
+
     - If the matrix is invertible, then it has an LU decomposition if and only
-    if all of its leading principal minors are non-zero (see
-    https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
-    leading principal minors of a matrix).
+      if all of its leading principal minors are non-zero (see
+      https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+      leading principal minors of a matrix).
     - If the matrix is singular (i.e., not invertible) and it has a rank of k
-    (i.e., it has k linearly independent columns), then it has an LU
-    decomposition if its first k leading principal minors are non-zero.
+      (i.e., it has k linearly independent columns), then it has an LU
+      decomposition if its first k leading principal minors are non-zero.
 
 This algorithm will simply attempt to perform LU decomposition on any square
 matrix and raise an error if no such decomposition exists.
@@ -25,6 +26,7 @@ def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray
     """
     Perform LU decomposition on a given matrix and raises an error if the matrix
     isn't square or if no such decomposition exists
+
     >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
     >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
     >>> lower_mat
@@ -45,7 +47,7 @@ def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray
     array([[ 4. ,  3. ],
            [ 0. , -1.5]])
 
-    # Matrix is not square
+    >>> # Matrix is not square
     >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
     >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
     Traceback (most recent call last):
@@ -54,14 +56,14 @@ def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray
     [[ 2 -2  1]
      [ 0  1  2]]
 
-    # Matrix is invertible, but its first leading principal minor is 0
+    >>> # Matrix is invertible, but its first leading principal minor is 0
     >>> matrix = np.array([[0, 1], [1, 0]])
     >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
     Traceback (most recent call last):
     ...
     ArithmeticError: No LU decomposition exists
 
-    # Matrix is singular, but its first leading principal minor is 1
+    >>> # Matrix is singular, but its first leading principal minor is 1
     >>> matrix = np.array([[1, 0], [1, 0]])
     >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
     >>> lower_mat
@@ -71,7 +73,7 @@ def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray
     array([[1., 0.],
            [0., 0.]])
 
-    # Matrix is singular, but its first leading principal minor is 0
+    >>> # Matrix is singular, but its first leading principal minor is 0
     >>> matrix = np.array([[0, 1], [0, 1]])
     >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
     Traceback (most recent call last):
diff --git a/linear_algebra/src/gaussian_elimination_pivoting.py b/linear_algebra/src/gaussian_elimination_pivoting.py
index efc1ddd64a2e..540f57b0cff6 100644
--- a/linear_algebra/src/gaussian_elimination_pivoting.py
+++ b/linear_algebra/src/gaussian_elimination_pivoting.py
@@ -6,17 +6,18 @@ def solve_linear_system(matrix: np.ndarray) -> np.ndarray:
     Solve a linear system of equations using Gaussian elimination with partial pivoting
 
     Args:
-    - matrix: Coefficient matrix with the last column representing the constants.
+      - `matrix`: Coefficient matrix with the last column representing the constants.
 
     Returns:
-    - Solution vector.
+      - Solution vector.
 
     Raises:
-    - ValueError: If the matrix is not correct (i.e., singular).
+      - ``ValueError``: If the matrix is not correct (i.e., singular).
 
     https://courses.engr.illinois.edu/cs357/su2013/lect.htm Lecture 7
 
     Example:
+
     >>> A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=float)
     >>> B = np.array([8, -11, -3], dtype=float)
     >>> solution = solve_linear_system(np.column_stack((A, B)))
diff --git a/linear_algebra/src/rank_of_matrix.py b/linear_algebra/src/rank_of_matrix.py
index 7ff3c1699a69..2c4fe2a8d1da 100644
--- a/linear_algebra/src/rank_of_matrix.py
+++ b/linear_algebra/src/rank_of_matrix.py
@@ -8,11 +8,15 @@
 def rank_of_matrix(matrix: list[list[int | float]]) -> int:
     """
     Finds the rank of a matrix.
+
     Args:
-        matrix: The matrix as a list of lists.
+        `matrix`: The matrix as a list of lists.
+
     Returns:
         The rank of the matrix.
+
     Example:
+
     >>> matrix1 = [[1, 2, 3],
     ...            [4, 5, 6],
     ...            [7, 8, 9]]
diff --git a/linear_algebra/src/schur_complement.py b/linear_algebra/src/schur_complement.py
index 7c79bb70abfc..74ac75e3fce2 100644
--- a/linear_algebra/src/schur_complement.py
+++ b/linear_algebra/src/schur_complement.py
@@ -12,13 +12,14 @@ def schur_complement(
 ) -> np.ndarray:
     """
     Schur complement of a symmetric matrix X given as a 2x2 block matrix
-    consisting of matrices A, B and C.
-    Matrix A must be quadratic and non-singular.
-    In case A is singular, a pseudo-inverse may be provided using
-    the pseudo_inv argument.
+    consisting of matrices `A`, `B` and `C`.
+    Matrix `A` must be quadratic and non-singular.
+    In case `A` is singular, a pseudo-inverse may be provided using
+    the `pseudo_inv` argument.
+
+    | Link to Wiki: https://en.wikipedia.org/wiki/Schur_complement
+    | See also Convex Optimization - Boyd and Vandenberghe, A.5.5
 
-    Link to Wiki: https://en.wikipedia.org/wiki/Schur_complement
-    See also Convex Optimization - Boyd and Vandenberghe, A.5.5
     >>> import numpy as np
     >>> a = np.array([[1, 2], [2, 1]])
     >>> b = np.array([[0, 3], [3, 0]])
diff --git a/linear_algebra/src/transformations_2d.py b/linear_algebra/src/transformations_2d.py
index b4185cd2848f..5dee59024752 100644
--- a/linear_algebra/src/transformations_2d.py
+++ b/linear_algebra/src/transformations_2d.py
@@ -3,13 +3,15 @@
 
 I have added the codes for reflection, projection, scaling and rotation 2D matrices.
 
+.. code-block:: python
+
     scaling(5) = [[5.0, 0.0], [0.0, 5.0]]
-  rotation(45) = [[0.5253219888177297, -0.8509035245341184],
-                  [0.8509035245341184, 0.5253219888177297]]
-projection(45) = [[0.27596319193541496, 0.446998331800279],
-                  [0.446998331800279, 0.7240368080645851]]
-reflection(45) = [[0.05064397763545947, 0.893996663600558],
-                  [0.893996663600558, 0.7018070490682369]]
+    rotation(45) = [[0.5253219888177297, -0.8509035245341184],
+                    [0.8509035245341184, 0.5253219888177297]]
+    projection(45) = [[0.27596319193541496, 0.446998331800279],
+                      [0.446998331800279, 0.7240368080645851]]
+    reflection(45) = [[0.05064397763545947, 0.893996663600558],
+                      [0.893996663600558, 0.7018070490682369]]
 """
 
 from math import cos, sin