diff --git a/src/scheme/ImplicitAcousticO2Solver.cpp b/src/scheme/ImplicitAcousticO2Solver.cpp
index 9351d911901caa8d75cc84d89f8fbb0a92fddd24..42af7235ffdd2cf122a2931b791605cd53e631f4 100644
--- a/src/scheme/ImplicitAcousticO2Solver.cpp
+++ b/src/scheme/ImplicitAcousticO2Solver.cpp
@@ -1707,15 +1707,44 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
     int number_iter = 0;
     do {
       number_iter++;
+      // first step
 
-      this->_computeGradJU_AU(u, p, pi, dt);
+      this->_computeGradJU_AU(u, p, pi, 0.5 * dt);
 
       NodeValue<const Rd> ur         = this->_getUr();
       NodeValuePerCell<const Rd> Fjr = this->_getFjr(ur);
+      // final step
 
-      // time n+1
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
 
+      CellValue<double> final_step_tau(m_mesh.connectivity());
+      CellValue<Rd> final_step_u = copy(u.cellValues());
+
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+          const auto& cell_nodes = cell_to_node_matrix[j];
+
+          Rd momentum_fluxes = zero;
+          for (size_t R = 0; R < cell_nodes.size(); ++R) {
+            const NodeId r = cell_nodes[R];
+            momentum_fluxes += Fjr(j, R);
+          }
+          const double dt_over_Mj = dt / m_Mj[j];
+          final_step_u[j] -= 0.5 * dt_over_Mj * momentum_fluxes;
+        });
+
+      for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+        double final_step_tau_fluxes = 0;
+        const auto& cell_nodes       = cell_to_node_matrix[j];
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          final_step_tau_fluxes += dot(m_Djr(j, R), ur[r]);
+        }
+        final_step_tau[j] = m_tau[j] + 0.5 * (dt / m_Mj[j]) * final_step_tau_fluxes;
+      };
+
+      // time n+1
+
       NodeValue<Rd> new_xr = copy(mesh->xr());
       parallel_for(
         mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
@@ -1725,7 +1754,7 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
       CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
 
       new_rho = copy(rho.cellValues());
-      new_u   = copy(u.cellValues());
+      new_u   = final_step_u;
       new_E   = copy(E.cellValues());
 
       parallel_for(
@@ -1740,8 +1769,8 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
             energy_fluxes += dot(Fjr(j, R), ur[r]);
           }
           const double dt_over_Mj = dt / m_Mj[j];
-          new_u[j] -= dt_over_Mj * momentum_fluxes;
-          new_E[j] -= dt_over_Mj * energy_fluxes;
+          new_u[j] -= 0.5 * dt_over_Mj * momentum_fluxes;
+          new_E[j] -= 0.5 * dt_over_Mj * energy_fluxes;
         });
 
       // update Djr
@@ -1786,7 +1815,7 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
           const NodeId r = cell_nodes[R];
           new_tau_fluxes += dot(m_Djr(j, R), ur[r]);
         }
-        new_tau[j] = m_tau[j] + (dt / m_Mj[j]) * new_tau_fluxes;
+        new_tau[j] = final_step_tau[j] + 0.5 * (dt / m_Mj[j]) * new_tau_fluxes;
         // std::cout << "new_tau(" << j << ")=" << new_tau[j] << '\n';
       }