diff --git a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl
index 61e202f6984..0f90ea22199 100644
--- a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl
+++ b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl
@@ -56,12 +56,12 @@ save_solution = SaveSolutionCallback(dt = 0.1,
 # Construct embedded order paired explicit Runge-Kutta method with 10 stages and 7 evaluation stages for given simulation setup.
 # Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used 
 # in calculating the polynomial coefficients in the ODE algorithm.
-ode_algorithm = Trixi.EmbeddedPairedRK3(10, 7, tspan, semi)
+ode_algorithm = Trixi.EmbeddedPairedRK3(10, 9, tspan, semi)
 
 # Calculate the CFL number for the given ODE algorithm and ODE problem (cfl_number calculate from dt_opt of the optimization of
 # b values in the Butcher tableau of the ODE algorithm).
 #cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
-stepsize_callback = StepsizeCallback() # Warisa: This number is quite small in contrast the other one from optimizing A
+stepsize_callback = StepsizeCallback()
 # I've tried using cfl of 1.5 and the error is very similar.
 
 callbacks = CallbackSet(summary_callback,
diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl
index f29c1cffee3..1f145fb76cb 100644
--- a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl
@@ -117,7 +117,7 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
         println("dt_opt_a = ", dt_opt_a)
         
 
-        error("b found.")
+        #error("b found.")
     end
     # Fill A-matrix in P-ERK style
     a_matrix = zeros(num_stage_evals - 2, 2)
@@ -370,7 +370,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator)
         end
 
         # Higher stages where the weight of b in the butcher tableau is zero
-        for stage in 2:(alg.num_stages - alg.num_stage_evals + 1)
+        for stage in 2:(alg.num_stages - alg.num_stage_evals)
             # Construct current state
             @threaded for i in eachindex(integrator.du)
                 integrator.u_tmp[i] = integrator.u[i] + alg.c[stage] * integrator.k1[i]
@@ -384,25 +384,28 @@ function step!(integrator::EmbeddedPairedRK3Integrator)
             end
         end
 
-        # k_(s-e+2)
+        # #k_(s-e+1) and k_(s-e+2)
         # Construct current state
-        @threaded for i in eachindex(integrator.du)
-            integrator.u_tmp[i] = integrator.u[i] +
-                                  alg.c[alg.num_stages - alg.num_stage_evals + 2] *
-                                  integrator.k1[i]
-        end
 
-        integrator.f(integrator.du, integrator.u_tmp, prob.p,
-                     integrator.t +
-                     alg.c[alg.num_stages - alg.num_stage_evals + 2] * integrator.dt)
+        for j in 1:2
+            @threaded for i in eachindex(integrator.du)
+                integrator.u_tmp[i] = integrator.u[i] +
+                                    alg.c[alg.num_stages - alg.num_stage_evals + 2] *
+                                    integrator.k1[i]
+            end
 
-        @threaded for i in eachindex(integrator.du)
-            integrator.k_higher[i] = integrator.du[i] * integrator.dt
-        end
+            integrator.f(integrator.du, integrator.u_tmp, prob.p,
+                        integrator.t +
+                        alg.c[alg.num_stages - alg.num_stage_evals + j] * integrator.dt)
 
-        @threaded for i in eachindex(integrator.u)
-            integrator.u[i] += integrator.k_higher[i] *
-                               alg.b[2]
+            @threaded for i in eachindex(integrator.du)
+                integrator.k_higher[i] = integrator.du[i] * integrator.dt
+            end
+
+            @threaded for i in eachindex(integrator.u)
+                integrator.u[i] += integrator.k_higher[i] *
+                                alg.b[j+1]
+            end
         end
 
         # Higher stages after num_stage_evals where b is non-zero
@@ -427,7 +430,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator)
 
             @threaded for i in eachindex(integrator.u)
                 integrator.u[i] += integrator.k_higher[i] *
-                                   alg.b[stage - alg.num_stages + alg.num_stage_evals]
+                                   alg.b[stage - alg.num_stages + alg.num_stage_evals+1]
             end
         end