diff --git a/ext/TrixiNLsolveExt.jl b/ext/TrixiNLsolveExt.jl
index 3db0e057f12..4516bdcf744 100644
--- a/ext/TrixiNLsolveExt.jl
+++ b/ext/TrixiNLsolveExt.jl
@@ -100,7 +100,7 @@ function Trixi.solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, num_stage_
         # Due to the nature of the nonlinear solver, different initial guesses can lead to 
         # small numerical differences in the solution.
 
-        x0 = convert(RealT, 0.1) .* rand(rng, RealT, num_stages - 2)
+        x0 = convert(RealT, 0.1) .* rand(rng, RealT, num_stage_evals - 2)
 
         sol = nlsolve(objective_function!, x0, method = :trust_region,
                       ftol = 4.0e-16, # Enforce objective up to machine precision
@@ -112,8 +112,8 @@ function Trixi.solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, num_stage_
         # and also c[i] - a[i, i-1] >= 0.0 since all coefficients should be non-negative
         # to avoid downwinding of numerical fluxes.
         is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown) &&
-                       all(!isnan(c[i] - a_unknown[i - 2]) &&
-                           c[i] - a_unknown[i - 2] >= 0 for i in eachindex(c) if i > 2)
+                       all(!isnan(c[i] - a_unknown[i - num_stage_evals]) &&
+                           c[i] - a_unknown[i - num_stage_evals] >= 0 for i in eachindex(c) if i > num_stage_evals)
 
         if is_sol_valid
             return a_unknown
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 3d63134aa6e..6540aced755 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
@@ -16,10 +16,10 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
     c = compute_c_coeffs(num_stages, cS2)
 
     # Initialize the array of our solution
-    a_unknown = zeros(num_stages - 2)
+    a_unknown = zeros(num_stage_evals - 2)
 
     # Special case of e = 3
-    if num_stages == 3
+    if num_stage_evals == 3
         a_unknown = [0.25] # Use classic SSPRK33 (Shu-Osher) Butcher Tableau
     else
         # Calculate coefficients of the stability polynomial in monomial form
@@ -45,8 +45,8 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
                                                     verbose)
     end
     # Fill A-matrix in P-ERK style
-    a_matrix = zeros(num_stages - 2, 2)
-    a_matrix[:, 1] = c[3:end]
+    a_matrix = zeros(num_stage_evals - 2, 2)
+    a_matrix[:, 1] = c[num_stages - num_stage_evals + 3:end] # issue
     a_matrix[:, 1] -= a_unknown
     a_matrix[:, 2] = a_unknown
 
@@ -303,14 +303,13 @@ function step!(integrator::EmbeddedPairedRK3Integrator)
         end
 
         # k_e (k at posoition 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]
+            integrator.u_tmp[i] = integrator.u[i] + alg.c[alg.num_stage_evals] * integrator.k1[i]
         end
 
         integrator.f(integrator.du, integrator.u_tmp, prob.p,
-                    integrator.t + alg.c[stage] * integrator.dt)
+                    integrator.t + alg.c[alg.num_stage_evals] * integrator.dt)
 
         @threaded for i in eachindex(integrator.du)
             integrator.k_higher[i] = integrator.du[i] * integrator.dt
@@ -326,9 +325,9 @@ function step!(integrator::EmbeddedPairedRK3Integrator)
             # Construct current state
             @threaded for i in eachindex(integrator.du)
                 integrator.u_tmp[i] = integrator.u[i] +
-                                      alg.a_matrix[stage - 2, 1] * #TODO: work on the indexing here and below
+                                      alg.a_matrix[stage - alg.num_stage_evals, 1]
                                       integrator.k1[i] +
-                                      alg.a_matrix[stage - 2, 2] *
+                                      alg.a_matrix[stage - alg.num_stage_evals, 2] *
                                       integrator.k_higher[i]
             end
 
@@ -341,24 +340,22 @@ 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 - 1]
+                                   alg.b[stage - alg.num_stage_evals + 1]
             end
         end
 
-        #TODO: Check if commenting this is equal to not commenting
-        #=
         # Last stage
         @threaded for i in eachindex(integrator.du)
             integrator.u_tmp[i] = integrator.u[i] +
-                                  alg.a_matrix[alg.num_stages - 2, 1] *
+                                  alg.a_matrix[alg.num_stages - alg.num_stage_evals - 1, 1] *
                                   integrator.k1[i] +
-                                  alg.a_matrix[alg.num_stages - 2, 2] *
+                                  alg.a_matrix[alg.num_stages - alg.num_stage_evals - 1, 2] *
                                   integrator.k_higher[i]
         end
 
         integrator.f(integrator.du, integrator.u_tmp, prob.p,
                      integrator.t + alg.c[alg.num_stages] * integrator.dt)
-        =#
+
         @threaded for i in eachindex(integrator.u)
             # add the contribution of the first stage
             integrator.u[i] += (1 - sum(alg.b)) * integrator.k1[i]