diff --git a/src/evolve.jl b/src/evolve.jl
index 854b01c2a5ef2c74abb7348515016dd7c3b6d3b0..c0675e5cd9acf6054f9543166dd707265e5f26f7 100644
--- a/src/evolve.jl
+++ b/src/evolve.jl
@@ -225,33 +225,24 @@ Stepper for Low Storage five-stage fourth-order Explicit Runge Kutta scheme.
 Taken from Hesthaven, Warburton (2007), section 3.4.
 """
 function ode_step!(u!, F, u0, t, dt, tmp_k, coeffs::COEFFS_LSERK4, callback_substep)
-
   @unpack nstages, a, b, c = coeffs
-
   kim1, ki = tmp_k
-  N = length(u!)
 
   for i = 1:nstages
-
     ti = t + c[i] * dt
     F(ki, i == 1 ? u0 : u!, ti)
-
-    @. ki *= dt
+    success = callback_substep(u!, ti, i+1) # +1 because u! is value for next stage
+    !success && return false
     if i == 1
-      @. u! = u0 + b[i] * ki
+      @. u! = u0 + b[i] * dt * ki
     else
       ai = a[i]
       @. ki += ai * kim1
-      @. u! += b[i] * ki
+      @. u! += b[i] * dt * ki
     end
-
-    success = callback_substep(u!, ti, i+1) # +1 because u! is value for next stage
-    !success && return false
-
     if i < nstages
       @. kim1 = ki
     end
-
   end # for
 
   return true
@@ -266,23 +257,23 @@ Stepper for Explicit Runge Kutta scheme.
 function ode_step!(u!, F, u0, t, dt, tmp_k, coeffs::TimeStepAlgorithm, callback_substep)
   @unpack nstages, a, b, c = coeffs
 
-  for n = 1:nstages
-    tn = t + c[n] * dt
+  for i = 1:nstages
+    ti = t + c[i] * dt
+    # assumes that a[1,:] are all zeros
+    F(tmp_k[i], i == 1 ? u0 : u!, ti)
+    i == nstages && break
+    success = callback_substep(u!, ti, i+1) # +1 because u! is value for next stage
+    !success && return false
     u! .= u0
-    for i = 1:n-1
-      u! .+= a[n, i] * tmp_k[i]
+    for j = 1:i
+      @. u! += a[i+1, j] * dt * tmp_k[j]
     end
-    F(tmp_k[n], u!, tn)
-    tmp_k[n] .*= dt
-    success = callback_substep(u!, tn, n)
-    !success && return false
   end
 
   u! .= u0
-  for n = 1:nstages
-    u! .+= b[n] * tmp_k[n]
+  for i = 1:nstages
+    @. u! += b[i] * dt * tmp_k[i]
   end
-
   success = callback_substep(u!, t+dt, nstages)
 
   return true
diff --git a/test/IntegrationTests/refs/grhd_bondi_infall_avmda/grhd_bondi_infall_avmda.toml b/test/IntegrationTests/refs/grhd_bondi_infall_avmda/grhd_bondi_infall_avmda.toml
new file mode 100644
index 0000000000000000000000000000000000000000..a5265ad9535316d2be6ab17e9c5e228524d5cb9b
--- /dev/null
+++ b/test/IntegrationTests/refs/grhd_bondi_infall_avmda/grhd_bondi_infall_avmda.toml
@@ -0,0 +1,33 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 0.02143872868864732
+polytrope_gamma = "$(4/3)"
+
+[GRHD]
+bc = "from_id"
+id = "bondi_accretion_infall"
+variables0d_analyze_error = ["ρ", "D"]
+variables1d_analyze_error = ["ρ", "D"]
+analyze_error_reference_solution = "bondi_accretion"
+c2p_enforce_causal_atm = false
+
+[Mesh]
+range = [ 1.8, 20.0 ]
+n = 7
+k = 16
+basis = "lgl"
+periodic = false
+
+[Output]
+variables1d = [ "D", "Sr", "Ï„" ]
+aligned_ts = "$(collect(range(1.0,10.0,step=1.0)))"
+
+[Evolution]
+cfl = 0.2
+tend = 10.0
+method = "rk4"
+
+[HRSC]
+method = "av"
+av_method = "mda"
+
diff --git a/test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5 b/test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5
new file mode 100644
index 0000000000000000000000000000000000000000..6f9227a415288f4c49e200f88648973114f66f4f
Binary files /dev/null and b/test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index cef00ea4e67e70d4b598be1657486bc0b1b5ff8b..8ae0756116439d53d34c62f4b10071f7eb9b301b 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -126,6 +126,9 @@ reltol1d = 1e-6
 variables0d  = [ "D_err_norm" ]
 variables1d  = [ "D_err", "D", "Sr", "τ", "p", "ρ", "ϵ", "vr" ]
 
+[grhd_bondi_infall_avmda]
+variables1d  = [ "D", "Sr", "Ï„", ]
+
 [fv_grhd_bondi_accretion]
 variables1d  = [ "D", "Sr", "τ", "p", "ρ", "ϵ", "vr" ]