diff --git a/src/one-body-sample.org b/src/one-body-sample.org
index 28b0996770676c06c0fe0eba26c89855956be7c6..bdb2f07bc3e44e8d9f7cd5aeff2c1355a9c29d0c 100644
--- a/src/one-body-sample.org
+++ b/src/one-body-sample.org
@@ -428,6 +428,7 @@
         double const time = tau * run;
         {
           stateUpdater->nextTimeStep();
+          timeSteppingScheme->nextTimeStep();
   
           VectorType ell(finestSize);
           assemble_neumann<GridType, GridView, SmallVector, P1Basis>
diff --git a/src/timestepping.org b/src/timestepping.org
index ed04802b99a272f3734b111d47c9ec37c8ab6bf2..a0d6f5bafc50e5118bb8ca515b281c77ade8c586 100644
--- a/src/timestepping.org
+++ b/src/timestepping.org
@@ -20,6 +20,7 @@
 #+end_src
 #+name: virtual_function_declaration
 #+begin_src c++
+  void virtual nextTimeStep();
   void virtual setup(VectorType const &, double, double, VectorType &,
                      VectorType &, MatrixType &);
   void virtual postProcess(VectorType const &);
@@ -57,6 +58,7 @@
   class TimeSteppingScheme
   {
   public:
+    void virtual nextTimeStep() = 0;
     void virtual
     setup(VectorType const &ell,
           double _tau,
@@ -69,6 +71,7 @@
     void virtual extractDisplacement(VectorType &displacement) const = 0;
     void virtual extractVelocity(VectorType &velocity) const = 0;
   };
+  
 #+end_src
 * TimeSteppingScheme: Implicit Euler
 #+name: ImplicitEuler.hh
@@ -116,6 +119,14 @@
       dirichletFunction(_dirichletFunction)
   {}
   
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
+  nextTimeStep()
+  {
+    ud_old = ud;
+    u_old = u;
+  }
+  
   template <class VectorType, class MatrixType, class FunctionType, int dim>
   void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
   setup(VectorType const &ell,
@@ -129,10 +140,6 @@
   
     tau = _tau;
   
-    // This function is only called once for every timestep
-    ud_old = ud;
-    u_old = u;
-  
     problem_rhs = ell;
     A.mmv(u_old, problem_rhs);
   
@@ -252,6 +259,15 @@
       dirichletFunction(_dirichletFunction)
   {}
   
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
+  nextTimeStep()
+  {
+    u_old_old = u_old;
+    ud_old = ud;
+    u_old = u;
+  }
+  
   template <class VectorType, class MatrixType, class FunctionType, int dim>
   void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
   setup(VectorType const &ell,
@@ -265,11 +281,6 @@
   
     tau = _tau;
   
-    // This function is only called once for every timestep
-    u_old_old = u_old;
-    ud_old = ud;
-    u_old = u;
-  
     switch (state) {
       // Perform an implicit Euler step since we lack information
     case NO_SETUP:
@@ -441,6 +452,15 @@
       dirichletFunction(_dirichletFunction)
   {}
   
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void Newmark<VectorType, MatrixType, FunctionType, dim>::
+  nextTimeStep()
+  {
+    udd_old = udd;
+    ud_old = ud;
+    u_old = u;
+  }
+  
   template <class VectorType, class MatrixType, class FunctionType, int dim>
   void Newmark<VectorType, MatrixType, FunctionType, dim>::
   setup(VectorType const &ell,
@@ -454,11 +474,6 @@
   
     tau = _tau;
   
-    // This function is only called once for every timestep
-    udd_old = udd;
-    ud_old = ud;
-    u_old = u;
-  
     problem_rhs = ell;
     B.usmv( 2.0/tau,  ud_old, problem_rhs);
     B.usmv(     1.0, udd_old, problem_rhs);