diff --git a/src/scheme/ImplicitAcousticO2Solver2steps.cpp b/src/scheme/ImplicitAcousticO2Solver2steps.cpp
index 6baec129bf385bb277eead2e0a03058371a05990..fb9cf391e84a7a74f3eeacf7562c33b1a3f41e97 100644
--- a/src/scheme/ImplicitAcousticO2Solver2steps.cpp
+++ b/src/scheme/ImplicitAcousticO2Solver2steps.cpp
@@ -1655,6 +1655,7 @@ class ImplicitAcousticO2SolverHandler2steps::ImplicitAcousticO2Solver final
     CellValue<Rd> new_u       = copy(u.cellValues());
     CellValue<double> new_E   = copy(E.cellValues());
     std::shared_ptr<const MeshType> new_mesh;
+    std::shared_ptr<const MeshType> mesh_1;
 
     m_Mj = [&]() {
       const CellValue<const double>& Vj = mesh_data.Vj();
@@ -1715,9 +1716,10 @@ class ImplicitAcousticO2SolverHandler2steps::ImplicitAcousticO2Solver final
       // second step
 
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
-      auto u_1                        = copy(u);
-      CellValue<double> tau_1(m_mesh.connectivity());
-      auto p_1 = copy(p);
+
+      auto u_inter = copy(u);
+      CellValue<double> tau_inter(m_mesh.connectivity());
+      auto p_inter = copy(p);
 
       parallel_for(
         mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
@@ -1729,7 +1731,7 @@ class ImplicitAcousticO2SolverHandler2steps::ImplicitAcousticO2Solver final
             momentum_fluxes += Fjr_1(j, R);
           }
           const double dt_over_Mj = dt / m_Mj[j];
-          u_1[j] -= 0.25 * dt_over_Mj * momentum_fluxes;
+          u_inter[j] -= 0.5 * dt_over_Mj * momentum_fluxes;
         });
 
       for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
@@ -1739,24 +1741,14 @@ class ImplicitAcousticO2SolverHandler2steps::ImplicitAcousticO2Solver final
           const NodeId r = cell_nodes[R];
           tau_fluxes_1 += dot(m_Djr(j, R), ur_1[r]);
         }
-        tau_1[j] = m_tau[j] + 0.25 * (dt / m_Mj[j]) * tau_fluxes_1;
-        p_1[j]   = std::pow((gamma[j] - 1) * exp(entropy[j] / Cv[j]), m_inv_gamma[j]) *
-                 pow(tau_1[j], gamma[j]);   // convertir tau en p par "p = (gamma-1)exp(S/C_v)tau^{-gamma}"
+        tau_inter[j] = m_tau[j] + 0.5 * (dt / m_Mj[j]) * tau_fluxes_1;
+        p_inter[j]   = std::pow((gamma[j] - 1) * exp(entropy[j] / Cv[j]), m_inv_gamma[j]) *
+                     pow(tau_inter[j], gamma[j]);   // convertir tau en p par "p = (gamma-1)exp(S/C_v)tau^{-gamma}"
       }
 
-      this->_computeGradJU_AU(u_1, p_1, pi, 0.25 * dt);
-
-      NodeValue<const Rd> ur_1_et_ur_2 = this->_getUr();
-
-      auto ur_2 = copy(ur_1_et_ur_2);
+      this->_computeGradJU_AU(u_inter, p_inter, pi, 0.25 * dt);
 
-      for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
-        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];
-          ur_2[r]        = ur_1_et_ur_2[r] - ur_1[r];
-        }
-      }
+      NodeValue<const Rd> ur_2 = this->_getUr();
 
       NodeValuePerCell<const Rd> Fjr_2 = this->_getFjr(ur_2);
 
@@ -1784,12 +1776,12 @@ class ImplicitAcousticO2SolverHandler2steps::ImplicitAcousticO2Solver final
           double energy_fluxes = 0;
           for (size_t R = 0; R < cell_nodes.size(); ++R) {
             const NodeId r = cell_nodes[R];
-            momentum_fluxes += Fjr_2(j, R);
-            energy_fluxes += dot(Fjr_2(j, R), ur_2[r]);
+            momentum_fluxes += Fjr_2(j, R) + Fjr_1(j, R);
+            energy_fluxes += dot(Fjr_2(j, R), ur_2[r]) + dot(Fjr_1(j, R), ur_1[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
@@ -1832,9 +1824,9 @@ class ImplicitAcousticO2SolverHandler2steps::ImplicitAcousticO2Solver final
         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];
-          new_tau_fluxes += dot(m_Djr(j, R), ur_2[r]);
+          new_tau_fluxes += dot(m_Djr(j, R), ur_2[r] + ur_1[r]);
         }
-        new_tau[j] = m_tau[j] + (dt / m_Mj[j]) * new_tau_fluxes;
+        new_tau[j] = m_tau[j] + 0.5 * (dt / m_Mj[j]) * new_tau_fluxes;
         // std::cout << "new_tau(" << j << ")=" << new_tau[j] << '\n';
       }