diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
new file mode 100644
index 0000000..29d872a
--- /dev/null
+++ b/.github/workflows/tests.yaml
@@ -0,0 +1,26 @@
+name: Run Tests
+
+on: [push]
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v2
+
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: 3.8
+
+    - name: Install dependencies
+      run: |
+        pip install -r requirements.txt
+      working-directory: ./UsingDorianFeatures
+
+    - name: Run tests
+      run: |
+        pytest
+      working-directory: ./UsingDorianFeatures
diff --git a/utils/lucas_heuristic.py b/utils/lucas_heuristic.py
new file mode 100644
index 0000000..2ac19ae
--- /dev/null
+++ b/utils/lucas_heuristic.py
@@ -0,0 +1,222 @@
+import sympy as sp
+import numpy as np
+from typing import List, Dict
+
+
+def newton_raphson(equations: List[sp.Expr],
+                   variables: List[sp.Symbol],
+                   initial_guess: List[float],
+                   tol: float = 1e-6,
+                   max_iter: int = 100,
+                   perturb_factor: float = 1e-5) -> List[float]:
+    """
+    Approximate the root of a set of multivariate polynomials
+    using the Newton-Raphson method.
+
+    Args:
+    - equations (list of sympy expressions): The system of equations to solve.
+    - variables (list of sympy symbols): The variables in the equations.
+    - initial_guess (list of float): Initial guess for the root.
+    - tol (float): Tolerance for convergence.
+    - max_iter (int): Maximum number of iterations.
+    - perturb_factor (float): Factor to perturb the root when the
+    Jacobian determinant is close to zero.
+
+    Returns:
+    - list of float: Approximated root.
+    """
+
+    # print(equations, variables, initial_guess)
+    m_vars = sp.Matrix(variables)
+    f = sp.Matrix(equations)
+    J = f.jacobian(m_vars)
+    root = sp.Matrix(initial_guess)
+
+    for i in range(max_iter):
+        f_eval = [float(expr.subs(list(zip(m_vars, root)))) for expr in f]
+
+        if all(abs(val) < tol for val in f_eval):
+            # print('max_iter', i)
+            return [float(val) for val in root]
+
+        J_eval = J.subs(list(zip(m_vars, root)))
+        determinant = J_eval.det()
+
+        # Set a counter
+        tries = 0
+        while abs(determinant) < tol and tries < max_iter:
+            tries += 1
+            # print(tries, max_iter)
+            # Perturb the root slightly
+            perturbation = np.random.rand(len(variables)) * perturb_factor
+            root += sp.Matrix(perturbation)
+            J_eval = J.subs(list(zip(m_vars, root)))
+            determinant = J_eval.det()
+
+        # Exit if the determinant is still to low
+        if tries >= max_iter:
+            break
+
+        delta_m_vars = -J_eval.inv() * sp.Matrix(f_eval)
+        root += delta_m_vars
+
+    raise ValueError("Newton-Raphson method did not converge within "
+                     "the specified maximum number of iterations.")
+
+
+def find_roots(equations: List[sp.Expr],
+               variables: List[sp.Symbol],
+               tol: float = 1e-6,
+               trials_factor: int = 5,
+               min_trials: int = 10) -> List[List[float]]:
+    """
+    Find roots of a set of multivariate polynomials
+    using the Newton-Raphson method
+    with automatic initial guess generation and
+    stopping criteria based on equation characteristics.
+
+    Args:
+    - equations (list of sympy expressions): The system of equations to solve.
+    - variables (list of sympy symbols): The variables in the equations.
+    - tol (float): Tolerance for convergence.
+    - trials_factor (int): The higher, the more times we look for new roots
+    - min_trials (int): Minimum number of times we look for new roots
+
+    Returns:
+    - list of lists of float: Approximated roots.
+    """
+
+    if len(equations) != len(variables):
+        raise ValueError("The number of independent equations should be equal"
+                         "to the number of variables to have finite roots.")
+
+    roots = []
+    attempts = 0
+
+    # Choose the range for the initial guesses
+    guess_range = 0
+    # Choose maximum number of attempts to find new roots
+    max_attempts = min_trials
+    for polynomial in equations:
+        # Compute the total degree of the polynomial
+        total_degree = polynomial.as_poly().total_degree()
+        if total_degree != 0:
+            abs_indep_coeff = abs(polynomial.as_coefficients_dict()[1])
+            guess_range = max(guess_range,
+                              2 * abs_indep_coeff ** (1 / total_degree))
+            max_attempts = total_degree * trials_factor
+
+    while attempts < max_attempts:
+        # Generate a random initial guess
+        initial_guess = [np.random.uniform(-guess_range, guess_range)
+                         for _ in variables]
+        # print('initial_guess', initial_guess)
+
+        try:
+            root = newton_raphson(equations, variables, initial_guess, tol)
+            # print('root', root)
+
+            # Check if the root is new (within tolerance)
+            is_new = all(
+                np.linalg.norm(np.array(root) - np.array(existing_root)) > tol
+                for existing_root in roots
+            )
+
+            if is_new:
+                roots.append(root)
+                attempts = 0  # Reset attempts if a new root is found
+            else:
+                # print('not new')
+                attempts += 1
+
+        except ValueError:
+            # Count the attempt if convergence error for this initial guess
+            attempts += 1
+
+    return roots
+
+
+def critical_points(polynomial: sp.Expr,
+                    trials_factor: int = 5) -> Dict[sp.Symbol, int]:
+    '''
+    Compute the critical points of the given polynomial with respect
+    to each of its variables.
+
+    Args:
+    - polynomial (sympy expression): The polynomial expression.
+    - trials_factor (int): The higher, the more times we look for new roots
+
+    Returns:
+    - dict: A dictionary containing the number of critical points
+            for each variable in the polynomial.
+    '''
+
+    # Find its variables
+    variables = list(polynomial.free_symbols)
+
+    # Store derivatives
+    derivatives = {}
+    for var in variables:
+        derivatives[var] = sp.diff(polynomial, var)
+
+    # Dictionary to store number of critical points
+    critical_points = dict()
+    for var in variables:
+        # print('var', var)
+        # List of equations containing:
+        # derivatives with respect to other variables
+        # and the original polynomial
+        equations = [derivatives[aux_var] for aux_var in variables
+                     if aux_var != var] + [polynomial]
+
+        # Use a root-finding function to find the critical points,
+        # that are the common roots of the equations
+        roots = find_roots(equations, variables, trials_factor=trials_factor)
+        # print('roots', roots)
+        critical_points[var] = len(roots)
+
+    return critical_points
+
+
+def count_critical_points(polynomials: List[sp.Expr]) -> Dict[sp.Symbol, int]:
+    """
+    Compute the total count of critical points
+    for each variable in a set of polynomials.
+
+    Args:
+    - polynomials (list of sympy expressions):
+    The set of multivariate polynomials to analyze.
+
+    Returns:
+    - dict: A dictionary where keys are variables
+    and values are the count of critical points.
+    """
+
+    # Find all unique variables present in the set of polynomials
+    unique_variables = set()
+    for polynomial in polynomials:
+        unique_variables.update(polynomial.free_symbols)
+
+    # Initialize a dictionary to store the total count of critical points
+    # for each variable
+    total_critical_points = {variable: 0 for variable in unique_variables}
+
+    # Iterate through each polynomial to find critical points
+    # and update the counts
+    for polynomial in polynomials:
+        for variable, count in critical_points(polynomial).items():
+            total_critical_points[variable] += count
+
+    return total_critical_points
+
+
+# Example usage:
+if __name__ == "__main__":
+    # Define a set of multivariate polynomials
+    x, y = sp.symbols('x y')
+    polynomials = [x**2 + y**2 - 4, x**3 * y - 4 * y**3]
+    polynomials = []
+
+    # Compute the total count of critical values for each variable
+    critical_values = count_critical_points(polynomials)
+    print(critical_values)  # Example output: {y: 3, x: 3}
diff --git a/utils/test_lucas_heuristic.py b/utils/test_lucas_heuristic.py
new file mode 100644
index 0000000..ede2360
--- /dev/null
+++ b/utils/test_lucas_heuristic.py
@@ -0,0 +1,184 @@
+import sympy as sp
+import unittest
+import time
+
+# Import the critical_points function from your module
+from lucas_heuristic import critical_points
+from lucas_heuristic import find_roots
+from lucas_heuristic import newton_raphson
+from lucas_heuristic import count_critical_points
+
+
+class TestCountCriticalPoints(unittest.TestCase):
+
+    def test_empty_input(self):
+        # Test with an empty list of polynomials
+        result = count_critical_points([])
+        expected = dict()
+        self.assertEqual(result, expected)
+
+    def test_single_polynomial(self):
+        # Test with a single polynomial with known critical points
+        x, y = sp.symbols('x y')
+        polynomial = x**2 + y**2 - 4
+        result = count_critical_points([polynomial])
+        expected = {x: 2, y: 2}  # Expected number of critical points
+        self.assertEqual(result, expected)
+
+    def test_multiple_polynomials(self):
+        # Test with multiple polynomials with known critical points
+        x, y = sp.symbols('x y')
+        polynomials = [x**2 + y**2 - 4,
+                       x**3 * y - 4 * y**3]
+        result = count_critical_points(polynomials)
+        expected = {x: 3, y: 3}  # Expected number of critical points
+        self.assertEqual(result, expected)
+
+
+class TestNewtonRaphson(unittest.TestCase):
+
+    def test_valid_input(self):
+        # Test with a valid set of equations, variables, and initial guess
+        x, y = sp.symbols('x y')
+        equations = [x**2 - 4,
+                     y**2 - 9]  # Two simple equations with roots (+-2, +-3)
+        variables = [x, y]
+        initial_guess = [1, 1]
+        root = newton_raphson(equations, variables, initial_guess)
+
+        # Ensure the root is one of the possibilities
+        self.assertAlmostEqual(root[0], 2, msg="The root given is not a root")
+        self.assertAlmostEqual(root[1], 3, msg="The root given is not a root")
+
+    def test_invalid_input(self):
+        # Test with invalid input where the number of equations
+        # and variables don't match
+        x, y, z = sp.symbols('x y z')
+        equations = [x**2 - 4, y**2 - 9]  # Two equations, but three variables
+        variables = [x, y, z]
+        initial_guess = [1, 1]
+
+        with self.assertRaises(ValueError):
+            newton_raphson(equations, variables, initial_guess)
+
+    def test_non_converging_case(self):
+        # Test a case where Newton-Raphson method
+        # doesn't converge within max_iter
+        x = sp.symbols('x')
+        equations = [x**2 + 1]  # No real roots
+        variables = [x]
+        initial_guess = [1]
+
+        with self.assertRaises(ValueError):
+            newton_raphson(equations, variables, initial_guess, max_iter=10)
+
+    def test_not_infinite_loop(self):
+        # It was observed that there can be an infinite loop
+        # this was fixed but example added here
+        # this test takes more than 10s
+        # if it enters the almost infinite loop
+
+        # Start measuring time
+        start_time = time.time()
+
+        x, y = sp.symbols('x y')
+        equations = [x**3, x**3*y - 4]
+        variables = [x, y]
+        initial_guess = [-1, 0]
+
+        with self.assertRaises(ValueError):
+            newton_raphson(equations, variables, initial_guess)
+
+        # Total time taken
+        execution_time = time.time() - start_time
+        # Assert that the execution time is less than 1 second
+        self.assertLess(execution_time, 1.0, "Test took longer than 1 second")
+
+
+class TestFindRoots(unittest.TestCase):
+
+    def test_valid_equations(self):
+        # Test with valid equations and variables
+        x, y = sp.symbols('x y')
+        equations = [x**2 + y**2 - 1,  # Circle equation
+                     x - y]            # diagonal line
+        variables = [x, y]
+        tol = 1e-6
+        roots = find_roots(equations, variables, tol=tol)
+
+        # Ensure the number of roots is correct (should be 2)
+        self.assertEqual(len(roots), 2)
+
+        # Ensure that they are the correct roots
+        sqrt2div2 = 0.7071068
+        for root in roots:
+            for coordinate in root:
+                self.assertAlmostEqual(abs(coordinate), sqrt2div2)
+
+    def test_invalid_input(self):
+        # Test with invalid input; number of equations
+        # and variables dont match
+        x, y, z = sp.symbols('x y z')
+        equations = [x**2 - 4, y**2 - 9]  # Two equations
+        variables = [x, y, z]             # but three variables
+
+        with self.assertRaises(ValueError):
+            find_roots(equations, variables)
+
+    def test_no_roots(self):
+        # Test with an invalid equation that cannot be solved
+        x, y = sp.symbols('x y')
+        equations = [x**2 + y**2 + 1,  # No real roots
+                     x - y]
+        variables = [x, y]
+        roots = find_roots(equations, variables)
+
+        # Ensure no roots are found (empty list)
+        self.assertEqual(len(roots), 0)
+
+
+class TestCriticalPoints(unittest.TestCase):
+
+    def test_single_variable(self):
+        # Test a polynomial with a single variable (e.g., x)
+        x = sp.symbols('x')
+        poly = x**2 + 3*x + 2
+        result = critical_points(poly)
+        expected = {x: 2}
+        self.assertDictEqual(result, expected)
+
+    def test_two_variables(self):
+        # Test a polynomial with two variables
+        x, y = sp.symbols('x y')
+        poly = x**2 + y**2 - 4
+        result = critical_points(poly)
+        expected = {x: 2, y: 2}
+        self.assertDictEqual(result, expected)
+
+    def test_no_independent_coefficient(self):
+        # Test a polynomial with no critical points
+        x, y = sp.symbols('x y')
+        poly = x**2 + y**2 * x
+        result = critical_points(poly)
+        expected = {x: 1, y: 1}
+        self.assertDictEqual(result, expected)
+
+    def test_no_critical_points(self):
+        # Test a polynomial with no critical points
+        x, y = sp.symbols('x y')
+        poly = x**2 + y**2 + 1
+        result = critical_points(poly)
+        expected = {x: 0, y: 0}
+        self.assertDictEqual(result, expected)
+
+    def test_complex_critical_points(self):
+        # Test a polynomial with complex critical points
+        x, y, z, w = sp.symbols('x y z w')
+        poly = x**4 + y**4 + z**4 + w**4 - 6*x**2 - 6*y**2 - 6*z**2 - 6*w**2
+        result = critical_points(poly, 10)
+        expected = {x: 1, y: 1, z: 1, w: 1}
+        self.assertDictEqual(result, expected)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/utils/trial.py b/utils/trial.py
new file mode 100644
index 0000000..22e7600
--- /dev/null
+++ b/utils/trial.py
@@ -0,0 +1,45 @@
+import sympy as sp
+
+
+def create_polynomial(monomials):
+    """
+    Create a polynomial from a list of monomials, where each monomial is represented as [degrees, coefficient].
+
+    Args:
+    - monomials (list of lists): List of monomials, where each monomial is represented as [degrees, coefficient].
+
+    Returns:
+    - sympy expression: The polynomial expression.
+    """
+
+    # Make sure the list of monomials is not empty
+    if monomials == []:
+        return 0
+
+    # Compute the number of variables
+    num_variables = len(monomials[0]) - 1
+
+    # Create symbolic variables in a loop and store them in a list
+    variables = [sp.symbols(f'x{i}') for i in range(num_variables)]
+
+    polynomial = sum(coefficient * sp.prod(var**degree
+                     for var, degree in zip(variables, degrees))
+                     for *degrees, coefficient in monomials)
+
+    return polynomial
+
+# Example usage:
+# if __name__ == "__main__":
+
+#     # List of monomials, each represented as [degrees, coefficient]
+#     monomials = [[1, 2, 0, 1], [1, 2, 3, 4]]
+
+#     # Create the polynomial
+#     polynomial = create_polynomial(monomials)
+
+#     # Print the polynomial
+#     print("Polynomial:", polynomial)
+
+dict1 = {'a': 2, 'b': 1, 'c': 3}
+dict2 = {'b': 1, 'a': 2, 'c': 3}
+assertDictEqual(dict1, dict2)
\ No newline at end of file