diff --git a/modmesh/app/euler1d.py b/modmesh/app/euler1d.py
index f1f79598e..3cc65b26b 100644
--- a/modmesh/app/euler1d.py
+++ b/modmesh/app/euler1d.py
@@ -530,7 +530,7 @@ def build_grid_figure(self):
 
         :return: FigureCanvas
         """
-        x = self.st.svr.coord[::2]
+        x = self.st.svr.coord[self.st.svr.xindices]
         fig = Figure()
         canvas = FigureCanvas(fig)
         ax = canvas.figure.subplots(3, 2)
@@ -574,7 +574,7 @@ def build_single_figure(self):
 
         :return: FigureCanvas
         """
-        x = self.st.svr.coord[::2]
+        x = self.st.svr.coord[self.st.svr.xindices]
         fig = Figure()
         canvas = FigureCanvas(fig)
         ax = canvas.figure.subplots()
@@ -743,24 +743,25 @@ def update_lines(self):
         :return: None
         """
         if self.use_grid_layout:
+            _s = self.st.svr.xindices
             self.density.update(adata=self.st.density_field,
-                                ndata=self.st.svr.density[::2])
+                                ndata=self.st.svr.density[_s])
             self.pressure.update(adata=self.st.pressure_field,
-                                 ndata=self.st.svr.pressure[::2])
+                                 ndata=self.st.svr.pressure[_s])
             self.velocity.update(adata=self.st.velocity_field,
-                                 ndata=self.st.svr.velocity[::2])
+                                 ndata=self.st.svr.velocity[_s])
             self.temperature.update(adata=self.st.temperature_field,
-                                    ndata=self.st.svr.temperature[::2])
+                                    ndata=self.st.svr.temperature[_s])
             self.internal_energy.update(adata=(self.st.internal_energy_field),
                                         ndata=(self.st.svr.
-                                               internal_energy[::2]))
+                                               internal_energy[_s]))
             self.entropy.update(adata=self.st.entropy_field,
-                                ndata=self.st.svr.entropy[::2])
+                                ndata=self.st.svr.entropy[_s])
         else:
             for name, is_selected, *_ in self.plot_config.state:
                 if is_selected:
                     eval(f'(self.{name}.update(adata=self.st.{name}_field,'
-                         f' ndata=self.st.svr.{name}[::2]))')
+                         f' ndata=self.st.svr.{name}[self.st.svr.xindices]))')
 
 
 class PlotArea(PuiInQt):
diff --git a/modmesh/onedim/euler1d.py b/modmesh/onedim/euler1d.py
index 90297d8a8..201deaf4d 100644
--- a/modmesh/onedim/euler1d.py
+++ b/modmesh/onedim/euler1d.py
@@ -24,10 +24,16 @@ class Euler1DSolver:
     method.
     """
 
-    def __init__(self, xmin, xmax, ncoord, time_increment=0.05):
+    def __init__(self, xmin, xmax, ncoord, time_increment=0.05,
+                 keep_edge=False):
         self._core = self.init_solver(xmin, xmax, ncoord, time_increment,
                                       gamma=1.4)
         # gamma is 1.4 for air.
+        _ = self.ncoord - 1
+        start = 0 if keep_edge else 2
+        stop = _ if keep_edge else (_ - 2)
+        num = (stop - start) // 2 + 1
+        self.xindices = np.linspace(start, stop, num, dtype='int32')
 
     def __getattr__(self, name):
         return getattr(self._core, name)
@@ -116,7 +122,7 @@ def __init__(self):
         self.svr = None
 
     def build_numerical(self, xmin, xmax, ncoord, time_increment=0.05,
-                        xdiaphragm=0.0):
+                        xdiaphragm=0.0, keep_edge=False):
         """
         After :py:meth:`build_constant` is done, optionally build the
         numerical solver :py:attr:`svr`.
@@ -133,7 +139,8 @@ def build_numerical(self, xmin, xmax, ncoord, time_increment=0.05,
 
         # Initialize the numerical solver.
         self.svr = Euler1DSolver(xmin, xmax, ncoord,
-                                 time_increment=time_increment)
+                                 time_increment=time_increment,
+                                 keep_edge=keep_edge)
 
         # Fill gamma.
         self.svr.gamma.fill(self.gamma)
@@ -310,8 +317,9 @@ def build_field(self, t, coord=None):
         :param coord: If None, take the coordinate from the numerical solver.
         :return: None
         """
-        if None is coord:
-            coord = self.svr.coord[::2]  # Use the numerical solver.
+        if coord is None:
+            # set the x-coordinate for numerical solver.
+            coord = self.svr.coord[self.svr.xindices]
         self.coord = coord.copy()  # Make a copy; no write back to argument.
 
         # Determine the zone location and the Boolean selection arrays.
diff --git a/tests/test_onedim_euler.py b/tests/test_onedim_euler.py
index 9fde27194..312598ebf 100644
--- a/tests/test_onedim_euler.py
+++ b/tests/test_onedim_euler.py
@@ -104,7 +104,7 @@ def test_field_without_numerical(self):
 
     def test_field_with_numerical(self):
         self.st.build_numerical(xmin=-1, xmax=1, ncoord=21,
-                                time_increment=0.05)
+                                time_increment=0.05, keep_edge=True)
         self.st.build_field(t=0.0)
         self.assertIsInstance(self.st.svr, euler1d.Euler1DSolver)
         self.assertEqual(len(self.st.coord), 11)