diff --git a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl
index 60f812b3e54..7b305bb2308 100644
--- a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl
+++ b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl
@@ -56,13 +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(16, 11, tspan, semi)
+ode_algorithm = Trixi.EmbeddedPairedRK3(16, 5, 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
-# I've tried using cfl of 1.5 and the error is very similar.
+stepsize_callback = StepsizeCallback()
 
 callbacks = CallbackSet(summary_callback,
                         alive_callback,
@@ -78,24 +77,3 @@ sol = Trixi.solve(ode, ode_algorithm,
 
 # Print the timer summary
 summary_callback()
-
-# Some function defined so that I can check if the second order condition is met. This will be removed later.
-function construct_b_vector(b_unknown, num_stages_embedded, num_stage_evals_embedded)
-    # Construct the b vector
-    b = [
-        b_unknown[1],
-        zeros(Float64, num_stages_embedded - num_stage_evals_embedded)...,
-        b_unknown[2:end]...,
-        0
-    ]
-    return b
-end
-
-b = construct_b_vector(ode_algorithm.b, ode_algorithm.num_stages - 1,
-                       ode_algorithm.num_stage_evals - 1)
-println("dot(b, c) = ", dot(b, ode_algorithm.c))
-println("sum(b) = ", sum(b))
-
-println("dt_opt_a = ", ode_algorithm.dt_opt_a)
-println("dt_opt_b = ", ode_algorithm.dt_opt_b)
-println("ratio = ", ode_algorithm.dt_opt_b / ode_algorithm.dt_opt_a * 100)
\ No newline at end of file
diff --git a/ext/TrixiConvexClarabelExt.jl b/ext/TrixiConvexClarabelExt.jl
deleted file mode 100644
index bedc05031f7..00000000000
--- a/ext/TrixiConvexClarabelExt.jl
+++ /dev/null
@@ -1,173 +0,0 @@
-# Package extension for adding Convex-based features to Trixi.jl
-module TrixiConvexClarabelExt
-
-# Required for coefficient optimization in P-ERK scheme integrators
-if isdefined(Base, :get_extension)
-    using Convex: MOI, solve!, Variable, minimize, evaluate, dot
-    using Clarabel: Optimizer
-else
-    # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl
-    using ..Convex: MOI, solve!, Variable, minimize, evaluate
-    using ..Clarabel: Optimizer
-end
-
-# Use functions that are to be extended and additional symbols that are not exported
-using Trixi: Trixi, @muladd
-
-# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
-# Since these FMAs can increase the performance of many numerical algorithms,
-# we need to opt-in explicitly.
-# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
-@muladd begin
-#! format: noindent
-
-# Compute new version of stability polynomial of the embedded scheme for paired explicit Runge-Kutta 
-# up to specified consistency order(p = 2), including contributions from free coefficients for higher
-# orders, and return the maximum absolute value
-function embedded_scheme_stability_polynomials!(pnoms,
-                                                num_stages_embedded,
-                                                num_stage_evals_embedded,
-                                                normalized_powered_eigvals_scaled,
-                                                a, b, c)
-
-    # Construct a full b coefficient vector #TODO: is there a way to not do this and just use b directly?
-    b_coeff = [
-        1 - sum(b),
-        zeros(Float64, num_stages_embedded - num_stage_evals_embedded)...,
-        b...,
-        0
-    ]
-    num_eig_vals = length(pnoms)
-
-    # Initialize with 1 + z
-    for i in 1:num_eig_vals
-        pnoms[i] = 1.0 + normalized_powered_eigvals_scaled[i, 1]
-    end
-    # z^2: b^T * c
-    #pnoms += dot(b, c) * normalized_powered_eigvals_scaled[:, 2]
-    pnoms += 0.5 * normalized_powered_eigvals_scaled[:, 2]
-
-    # Contribution from free coefficients
-    for i in 3:num_stage_evals_embedded
-        sum = 0.0
-        for j in (i + num_stages_embedded - num_stage_evals_embedded):num_stages_embedded
-            prod = 1.0
-            for k in (3 + j - i):j
-                prod *= a[k]
-            end
-            sum += prod * b_coeff[j] * c[j - i + 2]
-        end
-        pnoms += sum * normalized_powered_eigvals_scaled[:, i]
-    end
-
-    # For optimization only the maximum is relevant
-    return maximum(abs(pnoms))
-end
-
-function solve_b_butcher_coeffs_unknown_new(num_eig_vals, eig_vals,
-                                              num_stages, num_stage_evals,
-                                              num_stages_embedded,
-                                              num_stage_evals_embedded,
-                                              a_unknown, c, dtmax, dteps)
-    dtmin = 0.0
-    dt = -1.0
-    abs_p = -1.0
-
-    # Construct a full a coefficient vector
-    a = zeros(num_stages)
-    num_a_unknown = length(a_unknown)
-
-    for i in 1:num_a_unknown
-        a[num_stages - i + 1] = a_unknown[num_a_unknown - i + 1]
-    end
-
-    # Construct stability polynomial for each eigenvalue
-    pnoms = ones(Complex{Float64}, num_eig_vals, 1)
-
-    # There are e - 2 free variables for the stability polynomial of the embedded scheme
-    b = Variable(num_stage_evals - 2)
-
-    normalized_powered_eigvals = zeros(Complex{Float64}, num_eig_vals, num_stage_evals-1)
-
-    for j in 1:num_stage_evals-1
-        #fac_j = factorial(j)
-        for i in 1:num_eig_vals
-            #normalized_powered_eigvals[i, j] = eig_vals[i]^j / fac_j
-            # Try first without factorial normalization
-            normalized_powered_eigvals[i, j] = eig_vals[i]^j
-        end
-    end
-
-    normalized_powered_eigvals_scaled = similar(normalized_powered_eigvals)
-
-    # Bisection on timestep
-    while dtmax - dtmin > dteps
-        dt = 0.5 * (dtmax + dtmin)
-
-        # Compute stability polynomial for current timestep
-        for k in 1:num_stage_evals-1
-            dt_k = dt^k
-            for i in 1:num_eig_vals
-                normalized_powered_eigvals_scaled[i, k] = dt_k *
-                                                          normalized_powered_eigvals[i,
-                                                                                     k]
-            end
-        end
-
-        # Second-order constraint
-        # Since c[1] is always 0 we can ignore the contribution of b[1] and only account for the ones from other non-zero entries of b
-        constraints = [b >= 0,
-            2 * dot(b, c[(num_stages - num_stage_evals + 2):(num_stages - 1)]) == 1.0]
-
-        # Use last optimal values for b in (potentially) next iteration
-        problem = minimize(embedded_scheme_stability_polynomials!(pnoms,
-                                                                  num_stages_embedded,
-                                                                  num_stage_evals_embedded,
-                                                                  normalized_powered_eigvals_scaled,
-                                                                  a, b, c), constraints)
-
-        #=                                               
-        solve!(problem,
-               # Parameters taken from default values for EiCOS
-               MOI.OptimizerWithAttributes(Optimizer, "b" => 0.99,
-                                           "delta" => 2e-7,
-                                           "feastol" => 1e-9,
-                                           "abstol" => 1e-9,
-                                           "reltol" => 1e-9,
-                                           "feastol_inacc" => 1e-4,
-                                           "abstol_inacc" => 5e-5,
-                                           "reltol_inacc" => 5e-5,
-                                           "nitref" => 9,
-                                           "maxit" => 100,
-                                           "verbose" => 3); silent_solver = true)
-        =#
-
-        
-        solve!(problem,
-               MOI.OptimizerWithAttributes(Optimizer,
-                                           "tol_feas" => 1e-4);
-               silent_solver = true)
-        
-        print("dt: ", dt, "\n")
-        print("Problem value: ", problem.optval, "\n")
-        abs_p = problem.optval
-
-        if abs_p < 1
-            dtmin = dt
-        else
-            dtmax = dt
-        end
-    end
-
-    b_opt = evaluate(b)
-
-    # Catch case S = 3 (only one opt. variable)
-    if isa(b_opt, Number)
-        b_opt = [b_opt]
-    end
-
-    return b_opt, dt
-end
-end # @muladd
-
-end # module TrixiConvexClarabelExt
diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl
index c63a767f038..544438da871 100644
--- a/ext/TrixiConvexECOSExt.jl
+++ b/ext/TrixiConvexECOSExt.jl
@@ -3,11 +3,11 @@ module TrixiConvexECOSExt
 
 # Required for coefficient optimization in P-ERK scheme integrators
 if isdefined(Base, :get_extension)
-    using Convex: MOI, solve!, Variable, minimize, evaluate
+    using Convex: MOI, solve!, Variable, minimize, evaluate, sumsquares
     using ECOS: Optimizer
 else
     # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl
-    using ..Convex: MOI, solve!, Variable, minimize, evaluate
+    using ..Convex: MOI, solve!, Variable, minimize, evaluate, sumsquares
     using ..ECOS: Optimizer
 end
 
@@ -16,6 +16,7 @@ using LinearAlgebra: eigvals
 
 # Use functions that are to be extended and additional symbols that are not exported
 using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, solve_b_embedded, @muladd
+using Trixi: Trixi, undo_normalization!, bisect_stability_polynomial, compute_b_embedded_coeffs, @muladd
 
 # By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
 # Since these FMAs can increase the performance of many numerical algorithms,
@@ -294,6 +295,59 @@ function Trixi.solve_b_embedded(consistency_order, num_eig_vals, num_stage_evals
 
     return b_embedded, dt
 end
+
+function Trixi.compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c)
+
+    A = zeros(num_stage_evals-1, num_stage_evals)
+    rhs = [1, 1/2, embedded_monomial_coeffs...]
+
+    # Define the variables
+    b_embedded = Variable(num_stage_evals)  # the unknown coefficients we want to solve for
+    
+    # Initialize A matrix for the constraints
+    A[1, :] .= 1  # sum(b) = 1 constraint row
+    for i in 1:num_stage_evals-1
+        A[2, i+1] = c[num_stages - num_stage_evals + i]  # second order constraint: dot(b, c) = 1/2
+    end
+
+    # Fill the A matrix with other constraints based on monomial coefficients
+    for i in 3:(num_stage_evals-1)
+        for j in i+1:num_stage_evals
+            A[i, j] = c[num_stages - num_stage_evals + j - 2]
+            for k in 1:i-2
+                A[i, j] *= a_unknown[j - 1 - k]  # recursive dependence on `a_unknown`
+            end
+        end
+    end
+
+    println("A matrix of finding b")
+    display(A)
+
+    # Set up weights to prioritize the first two constraints
+    weights = [2, 2, ones(num_stage_evals-3)...]  # Heavier weight for the first two rows
+    weighted_residual = weights .* (A * b_embedded - rhs)  # Element-wise multiplication for weighting
+
+    # Set up the objective to minimize the weighted norm of the difference
+    problem = minimize(sumsquares(weighted_residual), [b_embedded >= 0])
+
+    solve!(problem, # Parameters taken from default values for EiCOS
+                                MOI.OptimizerWithAttributes(Optimizer, "gamma" => 0.99,
+                                "delta" => 2e-7,
+                                "feastol" => 1e-9,
+                                "abstol" => 1e-9,
+                                "reltol" => 1e-9,
+                                "feastol_inacc" => 1e-7,
+                                "abstol_inacc" => 5e-6,
+                                "reltol_inacc" => 5e-7,
+                                "nitref" => 9,
+                                "maxit" => 1000000,
+                                "verbose" => 3); silent_solver = true)
+    
+
+    ot = problem.optval
+    println("Optimal value of the objective function: ", ot)
+    return evaluate(b_embedded)
+end
 end # @muladd
 
 end # module TrixiConvexECOSExt
diff --git a/ext/TrixiNLsolveExt.jl b/ext/TrixiNLsolveExt.jl
index 319134e3700..ded4e7e281a 100644
--- a/ext/TrixiNLsolveExt.jl
+++ b/ext/TrixiNLsolveExt.jl
@@ -174,56 +174,6 @@ function Trixi.solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, num_stage_
 
     error("Maximum number of iterations ($max_iter) reached. Cannot find valid sets of coefficients.")
 end
-
-function Trixi.solve_b_butcher_coeffs_unknown!(b_unknown, num_stages, num_stage_evals, embedded_monomial_coeffs, c, a_unknown;
-    verbose, max_iter = 100000)
-
-    verbose = true
-
-    # Construct a full a coefficient vector
-    a = zeros(num_stages)
-    num_a_unknown = length(a_unknown)
-
-    for i in 1:num_a_unknown
-        a[num_stages - i + 1] = a_unknown[num_a_unknown - i + 1]
-    end
-
-    # Define the objective_function
-    function embedded_scheme_objective_function!(b_eq, x)
-        return EmbeddedPairedExplicitRK3_butcher_tableau_objective_function!(b_eq, x,
-                                                                     num_stages,
-                                                                     num_stage_evals,
-                                                                     embedded_monomial_coeffs,
-                                                                     c, a)
-    end
-
-    # RealT is determined as the type of the first element in monomial_coeffs to ensure type consistency
-    RealT = typeof(embedded_monomial_coeffs[1])
-
-    # To ensure consistency and reproducibility of results across runs, we use 
-    # a seeded random initial guess.
-    rng = StableRNG(55555)
-
-    # Due to the nature of the nonlinear solver, different initial guesses can lead to 
-    # small numerical differences in the solution.
-
-    for _ in 1:max_iter
-
-        # There is e-2 free variables of b of the embedded scheme
-        x0 = convert(RealT, 0.1) .* rand(rng, RealT, num_stage_evals - 2)
-
-        sol = nlsolve(embedded_scheme_objective_function!, x0, method = :trust_region,
-                    ftol = 4.0e-16, # Enforce objective up to machine precision
-                    iterations = 10^4, xtol = 1.0e-13, autodiff = :forward)
-
-        b_unknown = sol.zero # Retrieve solution (root = zero)
-
-        is_sol_valid = all(x -> !isnan(x) && x >= 0, b_unknown) && (sum(b_unknown) <= 1.0)
-
-        return b_unknown
-
-    end
-end
 end # @muladd
 
 end # module TrixiNLsolveExt
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 f1a59eb4459..6a592493d49 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
@@ -7,8 +7,57 @@ using DelimitedFiles: readdlm
 @muladd begin
 #! format: noindent
 
-# To be overwritten in TrixiConvexECOSExt.jl
-function solve_b_embedded end
+#=
+function compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c)
+
+    A = zeros(num_stage_evals - 1, num_stage_evals - 1)
+    b_embedded = zeros(num_stage_evals - 1)
+    rhs = [1, 1/2, embedded_monomial_coeffs...]
+
+    # sum(b) = 1
+    A[1, :] .= 1
+
+    # The second order constraint: dot(b,c) = 1/2
+    for i in 2:num_stage_evals - 1
+        A[2, i] = c[num_stages - num_stage_evals + i]
+    end
+
+    # Fill the A matrix
+    for i in 3:(num_stage_evals - 1)
+        # z^i
+        for j in i: (num_stage_evals - 1)
+            println("i = ", i, ", j = ", j)
+            println("[num_stages - num_stage_evals + j - 1] = ", num_stages - num_stage_evals + j - 1)
+            A[i,j] = c[num_stages - num_stage_evals + j - 1]
+            # number of times a_unknown should be multiplied in each power of z
+            for k in 1: i-2
+                # so we want to get from a[k] * ... i-2 times (1 time is already accounted for by c-coeff)
+                # j-1 - k + 1 = j - k
+                println("a_unknown at index: ", j - k)
+                A[i, j] *= a_unknown[j - k] # a[k] is in the same stage as b[k-1] -> since we also store b_1
+            end
+        end
+        #rhs[i] /= factorial(i)
+    end
+
+    display(A)
+
+    b_embedded = A \ rhs
+    return b_embedded
+end
+=#
+
+# Some function defined so that I can check if the second order condition is met. This will be removed later.
+function construct_b_vector(b_unknown, num_stages_embedded, num_stage_evals_embedded)
+    # Construct the b vector
+    b = [
+        b_unknown[1],
+        zeros(Float64, num_stages_embedded - num_stage_evals_embedded-1)...,
+        b_unknown[2:end]...,
+        0
+    ]
+    return b
+end
 
 # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3
 # using a list of eigenvalues
@@ -54,6 +103,12 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
 
         println("dt_opt_b = ", dt_opt_b)
 
+        b = compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c)
+        
+        b_full = construct_b_vector(b, num_stages - 1, num_stage_evals - 1)
+
+        println("dot(b, c) = ", dot(b_full, c))
+        println("sum(b) = ", sum(b_full))
         println("b: ", b)
         println("dt_opt_a = ", dt_opt_a)
         
@@ -311,7 +366,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]
@@ -325,25 +380,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
@@ -368,7 +426,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
 
diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
index 52d2d307db7..af396cb3ae7 100644
--- a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
@@ -16,5 +16,5 @@ include("polynomial_optimizer.jl")
 # such that hey can be exported from Trixi.jl and extended in the TrixiConvexECOSExt package
 # extension or by the NLsolve-specific code loaded by Requires.jl
 function solve_a_butcher_coeffs_unknown! end
-function solve_b_butcher_coeffs_unknown! end
+function compute_b_embedded_coeffs end
 end # @muladd