diff --git a/python_tsp/exact/brute_force.py b/python_tsp/exact/brute_force.py
index 5979bb8..04f0c6d 100644
--- a/python_tsp/exact/brute_force.py
+++ b/python_tsp/exact/brute_force.py
@@ -9,6 +9,7 @@
 
 def solve_tsp_brute_force(
     distance_matrix: np.ndarray,
+    starting_node: int = 0,
 ) -> Tuple[Optional[List], Any]:
     """Solve TSP to optimality with a brute force approach
 
@@ -18,6 +19,9 @@ def solve_tsp_brute_force(
         Distance matrix of shape (n x n) with the (i, j) entry indicating the
         distance from node i to j. It does not need to be symmetric
 
+    starting_node
+        Determines the starting node of the final permutation. Defaults to 0.
+
     Returns
     -------
     A permutation of nodes from 0 to n that produces the least total
@@ -33,14 +37,18 @@ def solve_tsp_brute_force(
     reducing the possibilities to (n - 1)!.
     """
 
-    # Exclude 0 from the range since it is fixed as starting point
-    points = range(1, distance_matrix.shape[0])
+    # Exclude `starting_node` from the range since it is fixed
+    other_points = [
+        node
+        for node in range(distance_matrix.shape[0])
+        if node != starting_node
+    ]
     best_distance = np.inf
     best_permutation = None
 
-    for partial_permutation in permutations(points):
+    for partial_permutation in permutations(other_points):
         # Remember to add the starting node before evaluating it
-        permutation = [0] + list(partial_permutation)
+        permutation = [starting_node] + list(partial_permutation)
         distance = compute_permutation_distance(distance_matrix, permutation)
 
         if distance < best_distance:
diff --git a/python_tsp/exact/dynamic_programming.py b/python_tsp/exact/dynamic_programming.py
index c8d8ffd..a465de4 100644
--- a/python_tsp/exact/dynamic_programming.py
+++ b/python_tsp/exact/dynamic_programming.py
@@ -7,6 +7,7 @@
 def solve_tsp_dynamic_programming(
     distance_matrix: np.ndarray,
     maxsize: Optional[int] = None,
+    starting_node: int = 0,
 ) -> Tuple[List, float]:
     """
     Solve TSP to optimality with dynamic programming
@@ -22,6 +23,9 @@ def solve_tsp_dynamic_programming(
         size for the recursion tree. Defaults to `None`, which essentially
         means "take as much space as needed".
 
+    starting_node
+        Determines the starting node of the final permutation. Defaults to 0.
+
     Returns
     -------
     permutation
@@ -89,16 +93,20 @@ def solve_tsp_dynamic_programming(
     ---------
     https://en.wikipedia.org/wiki/Held%E2%80%93Karp_algorithm#cite_note-5
     """
-    # Get initial set {1, 2, ..., tsp_size} as a frozenset because @lru_cache
-    # requires a hashable type
-    N = frozenset(range(1, distance_matrix.shape[0]))
+    # Get initial set {0, 1, 2, ..., tsp_size} \ {starting_node} as a frozenset
+    # because @lru_cache requires a hashable type
+    N = frozenset(
+        node
+        for node in range(distance_matrix.shape[0])
+        if node != starting_node
+    )
     memo: Dict[Tuple, int] = {}
 
     # Step 1: get minimum distance
     @lru_cache(maxsize=maxsize)
     def dist(ni: int, N: frozenset) -> float:
         if not N:
-            return distance_matrix[ni, 0]
+            return distance_matrix[ni, starting_node]
 
         # Store the costs in the form (nj, dist(nj, N))
         costs = [
@@ -110,11 +118,11 @@ def dist(ni: int, N: frozenset) -> float:
 
         return min_cost
 
-    best_distance = dist(0, N)
+    best_distance = dist(starting_node, N)
 
     # Step 2: get path with the minimum distance
-    ni = 0  # start at the origin
-    solution = [0]
+    ni = starting_node  # start at the origin
+    solution = [starting_node]
 
     while N:
         ni = memo[(ni, N)]
diff --git a/python_tsp/heuristics/lin_kernighan.py b/python_tsp/heuristics/lin_kernighan.py
index 592cbb0..d7dfc9e 100644
--- a/python_tsp/heuristics/lin_kernighan.py
+++ b/python_tsp/heuristics/lin_kernighan.py
@@ -143,6 +143,7 @@ def solve_tsp_lin_kernighan(
     x0: Optional[List[int]] = None,
     log_file: Optional[str] = None,
     verbose: bool = False,
+    starting_node: int = 0,
 ) -> Tuple[List[int], float]:
     """
     Solve the Traveling Salesperson Problem using the Lin-Kernighan algorithm.
@@ -163,6 +164,9 @@ def solve_tsp_lin_kernighan(
     verbose
         If true, prints algorithm status every iteration.
 
+    starting_node
+        Determines the starting node of the final permutation. Defaults to 0.
+
     Returns
     -------
     Tuple
@@ -174,7 +178,7 @@ def solve_tsp_lin_kernighan(
     Chapter 5, Section 5.3.2.1: Lin-Kernighan Neighborhood, Springer, 2023.
     """
     hamiltonian_cycle, hamiltonian_cycle_distance = setup_initial_solution(
-        distance_matrix=distance_matrix, x0=x0
+        distance_matrix=distance_matrix, x0=x0, starting_node=starting_node
     )
     num_vertices = distance_matrix.shape[0]
     vertices = list(range(num_vertices))
diff --git a/python_tsp/heuristics/local_search.py b/python_tsp/heuristics/local_search.py
index 75f261d..907b21a 100644
--- a/python_tsp/heuristics/local_search.py
+++ b/python_tsp/heuristics/local_search.py
@@ -21,6 +21,7 @@ def solve_tsp_local_search(
     max_processing_time: Optional[float] = None,
     log_file: Optional[str] = None,
     verbose: bool = False,
+    starting_node: int = 0,
 ) -> Tuple[List, float]:
     """Solve a TSP problem with a local search heuristic
 
@@ -47,6 +48,9 @@ def solve_tsp_local_search(
     verbose
         If true, prints algorithm status every iteration
 
+    starting_node
+        Determines the starting node of the final permutation. Defaults to 0.
+
     Returns
     -------
     A permutation of nodes from 0 to n - 1 that produces the least total
@@ -65,7 +69,9 @@ def solve_tsp_local_search(
         3. Repeat step 2 until all neighbors of `x` are tried and there is no
         improvement. Return `x`, `fx` as solution.
     """
-    x, fx = setup_initial_solution(distance_matrix, x0)
+    x, fx = setup_initial_solution(
+        distance_matrix, x0, starting_node=starting_node
+    )
     max_processing_time = max_processing_time or np.inf
 
     log_file_handler = (
diff --git a/python_tsp/heuristics/record_to_record.py b/python_tsp/heuristics/record_to_record.py
index 89a0a60..97be0f7 100644
--- a/python_tsp/heuristics/record_to_record.py
+++ b/python_tsp/heuristics/record_to_record.py
@@ -23,6 +23,7 @@ def solve_tsp_record_to_record(
     max_iterations: Optional[int] = None,
     log_file: Optional[str] = None,
     verbose: bool = False,
+    starting_node: int = 0,
 ):
     """
     Solve the traveling Salesperson Problem using a
@@ -48,6 +49,9 @@ def solve_tsp_record_to_record(
     verbose
         If true, prints algorithm status every iteration.
 
+    starting_node
+        Determines the starting node of the final permutation. Defaults to 0.
+
     Returns
     -------
     Tuple
@@ -60,7 +64,9 @@ def solve_tsp_record_to_record(
     """
     n = distance_matrix.shape[0]
     max_iterations = max_iterations or n
-    x, fx = setup_initial_solution(distance_matrix=distance_matrix, x0=x0)
+    x, fx = setup_initial_solution(
+        distance_matrix=distance_matrix, x0=x0, starting_node=starting_node
+    )
 
     log_file_handler = (
         open(log_file, "w", encoding="utf-8") if log_file else None
diff --git a/python_tsp/heuristics/simulated_annealing.py b/python_tsp/heuristics/simulated_annealing.py
index c19252e..a2c0c2d 100644
--- a/python_tsp/heuristics/simulated_annealing.py
+++ b/python_tsp/heuristics/simulated_annealing.py
@@ -24,6 +24,7 @@ def solve_tsp_simulated_annealing(
     max_processing_time: Optional[float] = None,
     log_file: Optional[str] = None,
     verbose: bool = False,
+    starting_node: int = 0,
 ) -> Tuple[List, float]:
     """Solve a TSP problem using a Simulated Annealing
     The approach used here is the one proposed in [1].
@@ -58,6 +59,9 @@ def solve_tsp_simulated_annealing(
     verbose {False}
         If true, prints algorithm status every iteration
 
+    starting_node
+        Determines the starting node of the final permutation. Defaults to 0.
+
     Returns
     -------
     A permutation of nodes from 0 to n - 1 that produces the least total
@@ -71,7 +75,9 @@ def solve_tsp_simulated_annealing(
     case studies. Springer Science & Business Media, 2006.
     """
 
-    x, fx = setup_initial_solution(distance_matrix, x0)
+    x, fx = setup_initial_solution(
+        distance_matrix, x0, starting_node=starting_node
+    )
     temp = _initial_temperature(distance_matrix, x, fx, perturbation_scheme)
     max_processing_time = max_processing_time or inf
     log_file_handler = (
diff --git a/python_tsp/utils/setup_initial_solution.py b/python_tsp/utils/setup_initial_solution.py
index 0e0fbc6..1cb9013 100644
--- a/python_tsp/utils/setup_initial_solution.py
+++ b/python_tsp/utils/setup_initial_solution.py
@@ -1,4 +1,4 @@
-from random import sample
+from random import sample, choice
 from typing import List, Optional, Tuple
 
 import numpy as np
@@ -6,8 +6,33 @@
 from .permutation_distance import compute_permutation_distance
 
 
+DEFAULT_STARTING_NODE = 0
+STARTING_NODE_OUTSIDE_BOUNDARIES_MSG = (
+    "Starting node outside limits [0, num_nodes]"
+)
+ENDING_NODE_OUTSIDE_BOUNDARIES_MSG = (
+    "Ending node outside limits [0, num_nodes]"
+)
+STARTING_ENDING_NODE_ARE_EQUAL_MSG = (
+    "Starting and ending nodes cannot be equal"
+)
+INVALID_SIZE_INITIAL_SOLUTION_MSG = (
+    "`x0` size does not match the number of nodes from distance matrix"
+)
+INVALID_INITIAL_SOLUTION_MSG = (
+    "`x0` does not contain a permutation of all nodes"
+)
+MISMATCH_INPUT_ARGUMENTS_MSG = (
+    "`x0` first and last node do not match with "
+    "`starting_node` or `ending_node`"
+)
+
+
 def setup_initial_solution(
-    distance_matrix: np.ndarray, x0: Optional[List] = None
+    distance_matrix: np.ndarray,
+    x0: Optional[List[int]] = None,
+    starting_node: Optional[int] = DEFAULT_STARTING_NODE,
+    ending_node: Optional[int] = None,
 ) -> Tuple[List[int], float]:
     """Return initial solution and its objective value
 
@@ -21,6 +46,14 @@ def setup_initial_solution(
         Permutation of nodes from 0 to n - 1 indicating the starting solution.
         If not provided, a random list is created.
 
+    starting_node
+        First node to appear in the permutation.
+
+    ending_node
+        Last node to appear in the permutation. Note, despite the name, the
+        TSP is still by definition a cycle, so we always go back from
+        `ending_node` to `starting_node`.
+
     Returns
     -------
     x0
@@ -30,10 +63,87 @@ def setup_initial_solution(
     fx0
         Objective value of x0
     """
+    _validate_input_arguments(distance_matrix, x0, starting_node, ending_node)
 
     if not x0:
-        n = distance_matrix.shape[0]  # number of nodes
-        x0 = [0] + sample(range(1, n), n - 1)  # ensure 0 is the first node
+        num_nodes = distance_matrix.shape[0]
+        x0 = _build_initial_permutation(
+            num_nodes, starting_node=starting_node, ending_node=ending_node
+        )
 
     fx0 = compute_permutation_distance(distance_matrix, x0)
     return x0, fx0
+
+
+def _validate_input_arguments(
+    distance_matrix: np.ndarray,
+    x0: Optional[List[int]],
+    starting_node: Optional[int],
+    ending_node: Optional[int],
+) -> None:
+    """Validate combination of input parameters"""
+    num_nodes = distance_matrix.shape[0]
+
+    if x0:
+        all_nodes = set(range(num_nodes))
+        if len(x0) != num_nodes:
+            raise ValueError(INVALID_SIZE_INITIAL_SOLUTION_MSG)
+        if set(x0) != all_nodes:
+            raise ValueError(INVALID_INITIAL_SOLUTION_MSG)
+        if starting_node is not None and x0[0] != starting_node:
+            raise ValueError(MISMATCH_INPUT_ARGUMENTS_MSG)
+        if ending_node is not None and x0[-1] != ending_node:
+            raise ValueError(MISMATCH_INPUT_ARGUMENTS_MSG)
+
+    if starting_node is not None:
+        if starting_node < 0 or starting_node >= num_nodes:
+            raise ValueError(STARTING_NODE_OUTSIDE_BOUNDARIES_MSG)
+
+    if ending_node is not None:
+        if ending_node < 0 or ending_node >= num_nodes:
+            raise ValueError(ENDING_NODE_OUTSIDE_BOUNDARIES_MSG)
+
+    if starting_node is not None and ending_node is not None:
+        if starting_node == ending_node:
+            raise ValueError(STARTING_ENDING_NODE_ARE_EQUAL_MSG)
+
+
+def _build_initial_permutation(
+    num_nodes: int,
+    starting_node: Optional[int] = DEFAULT_STARTING_NODE,
+    ending_node: Optional[int] = None,
+) -> List[int]:
+    """
+    Build a random list of integers from 0 to `num_nodes` - 1 guaranteeing the
+    initial node is `starting_node` and the last one is `ending_node`.
+
+    If not provided, `starting_node` is assumed as 0 and `ending_node` is taken
+    randomly.
+    """
+
+    all_nodes = list(range(num_nodes))
+    starting_node = (
+        starting_node if starting_node is not None else choice(all_nodes)
+    )
+
+    all_nodes_except_starting_node = [
+        node for node in range(num_nodes) if node != starting_node
+    ]
+    ending_node = (
+        ending_node
+        if ending_node is not None
+        else choice(all_nodes_except_starting_node)
+    )
+
+    all_nodes_without_extremes = [
+        node
+        for node in range(num_nodes)
+        if node != starting_node and node != ending_node
+    ]
+    x0 = (
+        [starting_node]
+        + sample(all_nodes_without_extremes, len(all_nodes_without_extremes))
+        + [ending_node]
+    )
+
+    return x0
diff --git a/tests/utils/test_setup_initial_solution.py b/tests/utils/test_setup_initial_solution.py
index e0f141e..4754ad1 100644
--- a/tests/utils/test_setup_initial_solution.py
+++ b/tests/utils/test_setup_initial_solution.py
@@ -1,6 +1,17 @@
 import pytest
 
+import numpy as np
+
 from python_tsp.utils import setup_initial_solution
+from python_tsp.utils.setup_initial_solution import (
+    DEFAULT_STARTING_NODE,
+    STARTING_NODE_OUTSIDE_BOUNDARIES_MSG,
+    ENDING_NODE_OUTSIDE_BOUNDARIES_MSG,
+    STARTING_ENDING_NODE_ARE_EQUAL_MSG,
+    INVALID_SIZE_INITIAL_SOLUTION_MSG,
+    INVALID_INITIAL_SOLUTION_MSG,
+    MISMATCH_INPUT_ARGUMENTS_MSG,
+)
 from tests.data import (
     distance_matrix1,
     distance_matrix2,
@@ -32,11 +43,157 @@ def test_setup_return_random_valid_solution(distance_matrix):
     """
     The setup outputs a random valid permutation if no initial solution
     is provided. This permutation must contain all nodes from 0 to n - 1
-    and start at 0 (the root).
+    and start at DEFAULT_STARTING_NODE.
     """
 
     x, fx = setup_initial_solution(distance_matrix)
 
     assert set(x) == set(range(distance_matrix.shape[0]))
-    assert x[0] == 0
+    assert x[0] == DEFAULT_STARTING_NODE
     assert fx
+
+
+def test_setup__random_starting_node_if_none():
+    """If explicitly set to `None`, `x0` can start at any node"""
+    num_nodes = 10
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    x0, fx0 = setup_initial_solution(
+        distance_matrix, starting_node=None
+    )
+
+    assert set(x0) == set(range(num_nodes))
+    assert fx0
+
+
+def test_setup__respects_starting_node():
+    num_nodes = 10
+    starting_node = 5
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    x0, fx0 = setup_initial_solution(
+        distance_matrix, starting_node=starting_node
+    )
+
+    assert set(x0) == set(range(num_nodes))
+    assert x0[0] == starting_node
+    assert fx0
+
+
+def test_setup__respects_ending_node():
+    num_nodes = 10
+    ending_node = 5
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    x0, fx0 = setup_initial_solution(distance_matrix, ending_node=ending_node)
+
+    assert set(x0) == set(range(num_nodes))
+    assert x0[0] == 0  # it is always 0 by default
+    assert x0[-1] == ending_node
+    assert fx0
+
+
+def test_setup__respects_starting_and_ending_node():
+    num_nodes = 10
+    starting_node = 7
+    ending_node = 5
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    x0, fx0 = setup_initial_solution(
+        distance_matrix, starting_node=starting_node, ending_node=ending_node
+    )
+
+    assert set(x0) == set(range(num_nodes))
+    assert x0[0] == starting_node
+    assert x0[-1] == ending_node
+    assert fx0
+
+
+def test_setup_input_validation__starting_node_outside_boundaries():
+    num_nodes = 10
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, starting_node=-1)
+    assert STARTING_NODE_OUTSIDE_BOUNDARIES_MSG in str(exc)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, starting_node=num_nodes)
+    assert STARTING_NODE_OUTSIDE_BOUNDARIES_MSG in str(exc)
+
+
+def test_setup_input_validation__ending_node_outside_boundaries():
+    num_nodes = 10
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, ending_node=-1)
+    assert ENDING_NODE_OUTSIDE_BOUNDARIES_MSG in str(exc)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, ending_node=num_nodes)
+    assert ENDING_NODE_OUTSIDE_BOUNDARIES_MSG in str(exc)
+
+
+def test_setup_input_validation__starting_ending_nodes_equal():
+    num_nodes = 10
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, starting_node=5, ending_node=5)
+    assert STARTING_ENDING_NODE_ARE_EQUAL_MSG in str(exc)
+
+
+def test_setup_input_validation__x0_has_wrong_size():
+    num_nodes = 10
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, x0=list(range(num_nodes + 1)))
+    assert INVALID_SIZE_INITIAL_SOLUTION_MSG in str(exc)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, x0=list(range(num_nodes - 1)))
+    assert INVALID_SIZE_INITIAL_SOLUTION_MSG in str(exc)
+
+
+def test_setup_input_validation__x0_does_not_contain_all_nodes():
+    num_nodes = 10
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+    x0 = list(range(num_nodes))
+    x0[5] = num_nodes + 5  # replace with a node outside 0:num_nodes-1
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, x0=x0, starting_node=0)
+    assert INVALID_INITIAL_SOLUTION_MSG in str(exc)
+
+
+def test_setup_input_validation__mismatch_arguments():
+    """
+    If `x0` and either `starting_node` or `ending_node` is provided, raise an
+    exception if they do not match.
+    """
+    num_nodes = 10
+    distance_matrix = np.random.rand(num_nodes, num_nodes)
+    x0 = list(range(num_nodes))
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, x0=x0, starting_node=5)
+    assert MISMATCH_INPUT_ARGUMENTS_MSG in str(exc)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(distance_matrix, x0=x0, ending_node=2)
+    assert MISMATCH_INPUT_ARGUMENTS_MSG in str(exc)
+
+    with pytest.raises(ValueError) as exc:
+        setup_initial_solution(
+            distance_matrix, x0=x0, starting_node=5, ending_node=0
+        )
+    assert MISMATCH_INPUT_ARGUMENTS_MSG in str(exc)
+
+    # The following ones pass
+    setup_initial_solution(distance_matrix, x0=x0, starting_node=0)
+    setup_initial_solution(distance_matrix, x0=x0, ending_node=num_nodes - 1)
+    setup_initial_solution(
+        distance_matrix, x0=x0, starting_node=0, ending_node=num_nodes - 1
+    )
diff --git a/tests/utils/test_starting_node.py b/tests/utils/test_starting_node.py
new file mode 100644
index 0000000..b0d3c8e
--- /dev/null
+++ b/tests/utils/test_starting_node.py
@@ -0,0 +1,58 @@
+import pytest
+
+from python_tsp.utils import setup_initial_solution
+from python_tsp.utils.setup_initial_solution import STARTING_NODE_TOO_LARGE_MSG
+from python_tsp.exact import (
+    solve_tsp_dynamic_programming,
+    solve_tsp_brute_force,
+    solve_tsp_branch_and_bound,
+)
+from python_tsp.heuristics import (
+    solve_tsp_local_search,
+    solve_tsp_simulated_annealing,
+    solve_tsp_lin_kernighan,
+    solve_tsp_record_to_record,
+)
+from tests.data import distance_matrix1
+
+
+@pytest.mark.parametrize("starting_node", [0, 1, 2, 3, 4])
+def test_initial_solution_respects_starting_node(starting_node):
+    x0, _ = setup_initial_solution(
+        distance_matrix1, starting_node=starting_node
+    )
+
+    # Ensure all nodes are contained in the solution and that
+    # the starting node is as requested
+    assert set(x0) == set(range(distance_matrix1.shape[0]))
+    assert x0[0] == starting_node
+
+
+def test_exception_is_raise_if_starting_node_is_too_large():
+    with pytest.raises(ValueError) as e:
+        x0, _ = setup_initial_solution(distance_matrix1, starting_node=999)
+
+    assert str(e.value) == STARTING_NODE_TOO_LARGE_MSG
+
+
+@pytest.mark.parametrize(
+    "solver",
+    [
+        solve_tsp_brute_force,
+        solve_tsp_dynamic_programming,
+        solve_tsp_local_search,
+        solve_tsp_simulated_annealing,
+        # These are not working yet
+        # solve_tsp_branch_and_bound,
+        # solve_tsp_lin_kernighan,
+        # solve_tsp_record_to_record,
+    ],
+)
+def test_solver_solution_respects_starting_node(solver):
+    starting_node = 3
+    x, _ = solver(
+        distance_matrix=distance_matrix1, starting_node=starting_node
+    )
+
+    assert set(x) == set(range(distance_matrix1.shape[0]))
+    assert x[0] == starting_node