diff --git a/src/bvpy/gbvp.py b/src/bvpy/gbvp.py
index 283756876fd552e57c156753bd31567c947930d0..5552e13960cb9c2d1ec89f1172144846937385d4 100644
--- a/src/bvpy/gbvp.py
+++ b/src/bvpy/gbvp.py
@@ -114,9 +114,12 @@ class GBVP(BVP):
         self.set_time(tmin)
         self.growth_scheme.dt = dt
 
+        current_solution = self.solution.copy(deepcopy=True)
+        current_growth = self.growth_scheme.get_growth().copy(deepcopy=True)
+
         if store_steps:
-            self.solution_steps[self.time.current] = self.solution.copy(deepcopy=True)
-            self.growth_steps[self.time.current] = self.growth_scheme.get_growth()
+            self.solution_steps[self.time.current] = current_solution
+            self.growth_steps[self.time.current] = current_growth
 
         if filename is not None:
             with fe.HDF5File(MPI_COMM_WORLD, filename, 'w') as file:
@@ -124,8 +127,8 @@ class GBVP(BVP):
                 #file.write(self.domain.cdata, "cell_labels")
                 #file.write(self.domain.bdata, "boundary_labels")
 
-                file.write(self.solution, f'solution_t{self.time.current}')
-                file.write(self.growth_steps[self.time.current], f'growth_t{self.time.current}')
+                file.write(current_solution, f'solution_t{self.time.current}')
+                file.write(current_growth, f'growth_t{self.time.current}')
 
             # Add parameters related to the domain : sub_domain / sub_boundary
             # if MPI_COMM_WORLD.rank == 0:  # h5py does not support parallal mode, use only the main process
@@ -146,14 +149,17 @@ class GBVP(BVP):
         while tmax - self.time.current >= dt:
             self.step(dt, **kwargs)
 
+            current_solution = self.solution.copy(deepcopy=True)
+            current_growth = self.growth_scheme.get_growth().copy(deepcopy=True)
+
             if store_steps:
-                self.solution_steps[self.time.current] = self.solution.copy(deepcopy=True)
-                self.growth_steps[self.time.current] = self.growth_scheme.get_growth().copy(deepcopy=True)
+                self.solution_steps[self.time.current] = current_solution
+                self.growth_steps[self.time.current] = current_growth
 
             if filename is not None:
                 with fe.HDF5File(MPI_COMM_WORLD, filename, 'a') as file:
-                    file.write(self.solution, f'solution_t{self.time.current}')
-                    file.write(self.growth_steps[self.time.current], f'growth_t{self.time.current}')
+                    file.write(current_solution, f'solution_t{self.time.current}')
+                    file.write(current_growth, f'growth_t{self.time.current}')
 
             self._progressbar.update(self.time.current)
             self._progressbar.show()
diff --git a/test/test_gbvp.py b/test/test_gbvp.py
index 6b430085160686749a41ddeff629762d8c055ea0..3361dadd602cf9a0efb57363a9c5ad8c0dc7c28b 100644
--- a/test/test_gbvp.py
+++ b/test/test_gbvp.py
@@ -35,7 +35,7 @@ class TestGBVP(unittest.TestCase):
     def test_run_with_file(self):
         gbvp = GBVP(self.domain, self.vform, self.bc, growth_scheme=self.growth_scheme)
         gbvp.info()
-        gbvp.run(0, 0.4, 0.2, filename=self.dir.name+'/solution.xdmf', linear_solver='mumps')
+        gbvp.run(0, 0.4, 0.2, filename=self.dir.name+'/solution.h5', linear_solver='mumps')
 
     def test_run_without_file(self):
         ibvp = GBVP(self.domain, self.vform, self.bc, growth_scheme=self.growth_scheme)