diff --git a/pina/geometry/cartesian.py b/pina/geometry/cartesian.py
index 6008589..ef1a6fe 100644
--- a/pina/geometry/cartesian.py
+++ b/pina/geometry/cartesian.py
@@ -36,7 +36,7 @@ class CartesianDomain(Location):
         :return: Spatial variables defined in ``__init__()``
         :rtype: list[str]
         """
-        return list(self.fixed_.keys()) + list(self.range_.keys())
+        return sorted(list(self.fixed_.keys()) + list(self.range_.keys()))
 
     def update(self, new_domain):
         """Adding new dimensions on the ``CartesianDomain``
diff --git a/pina/geometry/ellipsoid.py b/pina/geometry/ellipsoid.py
index 0764d81..e809bf6 100644
--- a/pina/geometry/ellipsoid.py
+++ b/pina/geometry/ellipsoid.py
@@ -78,7 +78,7 @@ class EllipsoidDomain(Location):
         :return: Spatial variables defined in '__init__()'
         :rtype: list[str]
         """
-        return list(self.fixed_.keys()) + list(self.range_.keys())
+        return sorted(list(self.fixed_.keys()) + list(self.range_.keys()))
 
     def is_inside(self, point, check_border=False):
         """Check if a point is inside the ellipsoid domain.
@@ -279,7 +279,7 @@ class EllipsoidDomain(Location):
             return _single_points_sample(n, variables)
 
         if variables == "all":
-            variables = list(self.range_.keys()) + list(self.fixed_.keys())
+            variables = self.variables
 
         if mode in ["random"]:
             return _Nd_sampler(n, mode, variables)
diff --git a/pina/geometry/simplex.py b/pina/geometry/simplex.py
index 2e78872..8a604b7 100644
--- a/pina/geometry/simplex.py
+++ b/pina/geometry/simplex.py
@@ -76,7 +76,7 @@ class SimplexDomain(Location):
 
     @property
     def variables(self):
-        return self._vertices_matrix.labels
+        return sorted(self._vertices_matrix.labels)
 
     def _build_cartesian(self, vertices):
         """
diff --git a/pina/problem/abstract_problem.py b/pina/problem/abstract_problem.py
index 6ebba6c..90a33d9 100644
--- a/pina/problem/abstract_problem.py
+++ b/pina/problem/abstract_problem.py
@@ -237,6 +237,9 @@ class AbstractProblem(metaclass=ABCMeta):
                 self.input_variables
             ):
                 self._have_sampled_points[location] = True
+                self.input_pts[location] = self.input_pts[location].extract(
+                    sorted(self.input_variables)
+                )
 
     def add_points(self, new_points):
         """
diff --git a/tests/test_problem.py b/tests/test_problem.py
index 8dcd499..f75fd58 100644
--- a/tests/test_problem.py
+++ b/tests/test_problem.py
@@ -67,12 +67,9 @@ class Poisson(SpatialProblem):
     truth_solution = poisson_sol
 
 
-# make the problem
-poisson_problem = Poisson()
-
-
 def test_discretise_domain():
     n = 10
+    poisson_problem = Poisson()
     boundaries = ['gamma1', 'gamma2', 'gamma3', 'gamma4']
     poisson_problem.discretise_domain(n, 'grid', locations=boundaries)
     for b in boundaries:
@@ -95,6 +92,7 @@ def test_discretise_domain():
 
 def test_sampling_few_variables():
     n = 10
+    poisson_problem = Poisson()
     poisson_problem.discretise_domain(n,
                                       'grid',
                                       locations=['D'],
@@ -103,20 +101,33 @@ def test_sampling_few_variables():
     assert poisson_problem._have_sampled_points['D'] is False
 
 
-# def test_sampling_all_args():
-#     n = 10
-#     poisson_problem.discretise_domain(n, 'grid', locations=['D'])
-
-# def test_sampling_all_kwargs():
-#     n = 10
-#     poisson_problem.discretise_domain(n=n, mode='latin', locations=['D'])
-
-# def test_sampling_dict():
-#     n = 10
-#     poisson_problem.discretise_domain(
-#         {'variables': ['x', 'y'], 'mode': 'grid', 'n': n}, locations=['D'])
-
-# def test_sampling_mixed_args_kwargs():
-#     n = 10
-#     with pytest.raises(ValueError):
-#         poisson_problem.discretise_domain(n, mode='latin', locations=['D'])
+def test_variables_correct_order_sampling():
+    n = 10
+    poisson_problem = Poisson()
+    poisson_problem.discretise_domain(n,
+                                      'grid',
+                                      locations=['D'],
+                                      variables=['x'])
+    poisson_problem.discretise_domain(n,
+                                      'grid',
+                                      locations=['D'],
+                                      variables=['y'])
+    assert poisson_problem.input_pts['D'].labels == sorted(
+        poisson_problem.input_variables)
+    
+    poisson_problem.discretise_domain(n,
+                                      'grid',
+                                      locations=['D'])
+    assert poisson_problem.input_pts['D'].labels == sorted(
+        poisson_problem.input_variables)
+    
+    poisson_problem.discretise_domain(n,
+                                      'grid',
+                                      locations=['D'],
+                                      variables=['y'])
+    poisson_problem.discretise_domain(n,
+                                      'grid',
+                                      locations=['D'],
+                                      variables=['x'])
+    assert poisson_problem.input_pts['D'].labels == sorted(
+        poisson_problem.input_variables)
\ No newline at end of file