diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
index b11057c77c..ba4c51c0d0 100644
--- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -166,17 +166,18 @@ struct AdaptiveRadau{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
  new_W_γdt_cutoff::C2
  controller::Symbol
  step_limiter!::StepLimiter
+ num_stages::Int
 end
 
 function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      standardtag = Val{true}(), concrete_jac = nothing,
-     diff_type = Val{:forward},
+     diff_type = Val{:forward}, num_stages = 5, 
      linsolve = nothing, precs = DEFAULT_PRECS,
      extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
      new_W_γdt_cutoff = 1 // 5,
      controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
      step_limiter! = trivial_limiter!)
- RadauIIA9{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+ AdaptiveRadau{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
      typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
      typeof(κ), typeof(fast_convergence_cutoff),
      typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
@@ -188,6 +189,6 @@ function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      fast_convergence_cutoff,
      new_W_γdt_cutoff,
      controller,
-     step_limiter!)
+     step_limiter!, num_stages)
 end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index cf8bcd5c5b..40f6631765 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -468,7 +468,7 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits},
         Convergence, alg.step_limiter!)
 end
 
-mutable struct adaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType, S} <:
+mutable struct adaptiveRadauConstantCache{S, F, Tab, Tol, Dt, U, JType} <:
     OrdinaryDiffEqConstantCache
 uf::F
 tab::Tab
@@ -482,16 +482,16 @@ status::NLStatus
 J::JType
 end
 
-function alg_cache(alg::adaptiveRadau, s :: Int64, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
     ::Type{uBottomEltypeNoUnits},
     ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
     ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
 uf = UDerivativeWrapper(f, t, p)
 uToltype = constvalue(uBottomEltypeNoUnits)
-tab = adaptiveRadau(uToltype, constvalue(tTypeNoUnits), s)
+tab = adaptiveRadau(uToltype, constvalue(tTypeNoUnits), alg.num_stages)
 
-cont = Vector{typeof(u)}(undef, s-1)
-for i in 1:s-1
+cont = Vector{typeof(u)}(undef, num_stages - 1)
+for i in 1: (num_stages - 1)
     cont[i] = zero(u)
 end
 
@@ -539,7 +539,7 @@ mutable struct adaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, JType,
     step_limiter!::StepLimiter
 end
 
-function alg_cache(alg::adaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
     ::Type{uBottomEltypeNoUnits},
     ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
     ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
@@ -549,34 +549,34 @@ tab = RadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits))
 
 κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
 
-z = Vector{typeof(u)}(undef, s)
-w = Vector{typeof(u)}(undef, s)
+z = Vector{typeof(u)}(undef, num_stages)
+w = Vector{typeof(u)}(undef, num_stages)
 for i in 1:s
     z[i] = w[i] = zero(u)
 end
 
 dw1 = zero(u)
 ubuff = zero(u)
-dw2 = Vector{typeof(u)}(undef, floor(Int, s/2))
-for i in 1 : floor(Int, s/2)
+dw2 = Vector{typeof(u)}(undef, floor(Int, num_stages/2))
+for i in 1 : floor(Int, num_stages/2)
     dw2[i] = similar(u, Complex{eltype(u)})
     recursivefill!(dw[i], false)
 end
-cubuff = Vector{typeof(u)}(undef, floor(Int, s/2))
-for i in 1 :floor(Int, s/2)
+cubuff = Vector{typeof(u)}(undef, floor(Int, num_stages/2))
+for i in 1 :floor(Int, num_stages/2)
     cubuff[i] = similar(u, Complex{eltype(u)})
     recursivefill!(cubuff[i], false)
 end
 
-cont = Vector{typeof(u)}(undef, s-1)
-for i in 1:s-1
+cont = Vector{typeof(u)}(undef, num_stages - 1)
+for i in 1: (num_stages - 1)
     cont[i] = zero(u)
 end
 
 fsalfirst = zero(rate_prototype)
-fw = Vector{typeof(rate_prototype)}(undef, s)
-k = Vector{typeof(rate_prototype)}(undef, s)
-for i in 1:s
+fw = Vector{typeof(rate_prototype)}(undef, num_stages)
+k = Vector{typeof(rate_prototype)}(undef, num_stages)
+for i in 1: num_stages
     k[i] = fw[i] = zero(rate_prototype)
 end
 
@@ -584,16 +584,16 @@ J, W1 = build_J_W(alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val(true))
 if J isa AbstractSciMLOperator
     error("Non-concrete Jacobian not yet supported by RadauIIA5.")
 end
-W2 = vector{typeof(Complex{W1})}(undef, floor(Int, s/2))
-for i in 1 : floor(Int, s/2)
+W2 = vector{typeof(Complex{W1})}(undef, floor(Int, num_stages/2))
+for i in 1 : floor(Int, num_stages/2)
     W2[i] = similar(J, Complex{eltype(W1)})
     recursivefill!(w2[i], false)
 end
 
 du1 = zero(rate_prototype)
 
-tmp = Vector{typeof(u)}(undef, binomial(s,2))
-for i in 1 : binomial(s,2)
+tmp = Vector{typeof(u)}(undef, binomial(num_stages,2))
+for i in 1 : binomial(num_stages,2)
     tmp[i] = zero(u)
 end
 
@@ -606,8 +606,8 @@ linprob = LinearProblem(W1, _vec(ubuff); u0 = _vec(dw1))
 linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
     assumptions = LinearSolve.OperatorAssumptions(true))
 
-linsolve2 = Vector{typeof(linsolve1)}(undef, floor(Int, s/2))
-for i in 1 : floor(int, s/2)
+linsolve2 = Vector{typeof(linsolve1)}(undef, floor(Int, num_stages/2))
+for i in 1 : floor(int, num_stages/2)
     linprob = LinearProblem(W2[i], _vec(cubuff[i]); u0 = _vec(dw2[i]))
     linsolve2 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
     assumptions = LinearSolve.OperatorAssumptions(true))
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
index afd2392a48..51313203a6 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
@@ -1341,9 +1341,9 @@ end
 end
 
 @muladd function perform_step!(integrator, cache::adaptiveRadauConstantCache,
-    repeat_step = false, s::Int64)
+    repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack T, TI, γ, α, β, c, e= cache.tab 
+    @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab 
     @unpack κ, cont = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
@@ -1357,39 +1357,39 @@ end
     γdt, αdt, βdt = γ / dt, α ./ dt, β ./ dt
 
     J = calc_J(integrator, cache)
-    LU = Vector{Any}(undef, (s + 1) / 2)
+    LU = Vector{Any}(undef, (num_stages + 1) / 2)
     if u isa Number
         LU[1] = -γdt * mass_matrix + J
-        for i in 2 : (s + 1) / 2
+        for i in 2 : (num_stages + 1) / 2
             LU[i] = -(α[i - 1]dt + β[i - 1]dt * im) * mass_matrix + J
         end
     else
         LU[1] = lu(-γdt * mass_matrix + J)
-        for i in 2 : (s + 1) / 2
+        for i in 2 : (num_stages + 1) / 2
             LU[i] = lu(-(α[i - 1]dt + β[i - 1]dt * im) * mass_matrix + J)
         end
     end
     integrator.stats.nw += 1
-
+    z = w = Vector{BigFloat}(undef, num_stages)
     if integrator.iter == 1 || integrator.u_modified || alg.extrapolant == :constant
         cache.dtprev = one(cache.dtprev)
-        for i in 1:s
+        for i in 1 : num_stages
             z[i] = w[i] = map(zero, u)
         end
-        for i in 1:(s-1)
+        for i in 1 : (num_stages - 1)
             cont[i] = map(zero, u)
         end
     else 
-        c' = Vector{eltype(u)}(undef, s) #time stepping
-        c'[s] = dt / cache.dtprev
-        for i in 1 : s-1
-            c'[i] = c[i] * c'[s]
+        c' = Vector{eltype(u)}(undef, num_stages) #time stepping
+        c'[num_stages] = dt / cache.dtprev
+        for i in 1 : num_stages - 1
+            c'[i] = c[i] * c'[num_stages]
         end
-        for i in 1 : s # collocation polynomial
-            z[i] = @.. cont[s-1] * (c'[i] - c[1] + 1) + cont[s-1]
+        for i in 1 : num_stages # collocation polynomial
+            z[i] = @.. cont[num_stages-1] * (c'[i] - c[1] + 1) + cont[num_stages - 1]
             j = s - 2
             while j > 0
-                z[i] = @.. z[i] * (c'[i] - c[s-j] + 1) + cont[j]
+                z[i] = @.. z[i] * (c'[i] - c[num_stages-j] + 1) + cont[j]
             end
             z[i] = @.. z[i] * c'[i]
         end
@@ -1406,49 +1406,49 @@ end
         integrator.stats.nnonliniter += 1
 
         # evaluate function
-        ff = Vector{eltype(u)}(undef, s)
-        for i in 1:s
+        ff = Vector{eltype(u)}(undef, num_stages)
+        for i in 1 : num_stages
             ff[i] = f(uprev + z[i], p, t + c[i] * dt)
         end    
         integrator.stats.nf += 5
 
         fw = @.. TI * ff
-        Mw = Vector{eltype(u)}(undef, s)
+        Mw = Vector{eltype(u)}(undef, num_stages)
         if mass_matrix isa UniformScaling # `UniformScaling` doesn't play nicely with broadcast
-            for i in 1:s
-                Mw[i] = @.. mass_matrix.λ * w[i] #scaling by eigenvalues
+            for i in 1 : num_stages
+                Mw[i] = @.. mass_matrix.λ * w[i] #scaling by eigenvalue
             end
         else
             Mw = mass_matrix * w #standard multiplication
         end
 
-        rhs = Vector{eltype(u)}(undef, s)
-        rhs[1] = @.. fw[1]-γdt * Mw[1]
+        rhs = Vector{eltype(u)}(undef, num_stages)
+        rhs[1] = @.. fw[1] - γdt * Mw[1]
         i = 2
-        while i <= s #block by block multiplication
-            rhs[i] = @.. fw[i] - α[i/2]dt * Mw[i] + β[i/2]dt * Mw[i + 1]
-            rhs[i + 1] = @.. fw[i + 1] - β[i/2]dt * Mw[i] - α[i/2]dt * Mw[i + 1]
+        while i <= num_stages #block by block multiplication
+            rhs[i] = @.. fw[i] - αdt[i/2] * Mw[i] + βdt[i/2] * Mw[i + 1]
+            rhs[i + 1] = @.. fw[i + 1] - βdt[i/2] * Mw[i] - αdt[i/2] * Mw[i + 1]
             i += 2
         end 
 
-        dw = Vector{eltype(u)}(undef, s)
+        dw = Vector{eltype(u)}(undef, num_stages)
         dw[1] = LU1 \ rhs[1]
-        for i in 2 : (s + 1) / 2 
+        for i in 2 : (num_stages + 1) / 2 
             tmp = LU[i] \ (@.. rhs[2 * i - 2] + rhs[2 * i - 1] * im)
             dw[2 * i - 2] = real(tmp)
             dw[2 * i - 1] = imag(tmp)
         end
-        integrator.stats.nlsolve += (s + 1) / 2
+        integrator.stats.nlsolve += (num_stages + 1) / 2
 
         # compute norm of residuals
         iter > 1 && (ndwprev = ndw)
-        atmp = Vector{eltype(u)}(undef, s)
-        for i in 1:s
+        atmp = Vector{eltype(u)}(undef, num_stages)
+        for i in 1 : num_stages
             atmp[i] = calculate_residuals(dw[i], uprev, u, atol, rtol, internalnorm, t)
         end
 
         ndw = 0
-        for i in 1:s
+        for i in 1 : num_stages
             ndw = ndw + internalnorm(atmp[i], t)
         end
         # check divergence (not in initial step)
@@ -1487,8 +1487,8 @@ end
     cache.ηold = η
     cache.iter = iter
 
-    u = @.. uprev + z[s]
-
+    u = @.. uprev + z[num_stages]
+    #=
     if adaptive
         edt = e ./ dt
         tmp = @.. dot(edt, z)
@@ -1508,13 +1508,13 @@ end
             integrator.EEst = internalnorm(atmp, t)
         end
     end
-
+    =#
     if integrator.EEst <= oneunit(integrator.EEst)
         cache.dtprev = dt
         if alg.extrapolant != :constant
-            derivatives = Matrix{eltype(u)}(undef, s-1, s-1)
-            for i in 1 : (s - 1)
-                for j in i : (s-1)
+            derivatives = Matrix{eltype(u)}(undef, num_stages - 1, num_stages - 1)
+            for i in 1 : (num_stages - 1)
+                for j in i : (num_stages - 1)
                     if i == 1
                         derivatives[i, j] = @.. (z[i] - z[i + 1]) / (c[i] - c[i + 1]) #first derivatives
                     else
@@ -1522,8 +1522,8 @@ end
                     end
                 end
             end
-            for i in 1 : (s-1)
-                cache.cont[i] = derivatives[i, s - 1]
+            for i in 1 : (num_stages - 1)
+                cache.cont[i] = derivatives[i, num_stages - 1]
             end
         end
     end
@@ -1538,7 +1538,7 @@ end
 
 @muladd function perform_step!(integrator, cache::adaptiveRadauCache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p, fsallast, fsalfirst = integrator
-    @unpack T, TI, γ, α, β, c, e= cache.tab 
+    @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab 
     @unpack κ, cont, z, w = cache
     @unpack dw1, ubuff, dw2, cubuff1, cubuff2 = cache
     @unpack k, fw, J, W1, W2 = cache
@@ -1556,7 +1556,7 @@ end
     if (new_W = do_newW(integrator, alg, new_jac, cache.W_γdt))
         @inbounds for II in CartesianIndices(J)
             W1[II] = -γdt * mass_matrix[Tuple(II)...] + J[II]
-            for i in 1 : (s - 1) / 2
+            for i in 1 : (num_stages - 1) / 2
                 W2[i][II] = -(α[i]dt + β[i]dt * im) * mass_matrix[Tuple(II)...] + J[II]
             end        
         end
@@ -1567,23 +1567,23 @@ end
     if integrator.iter == 1 || integrator.u_modified || alg.extrapolant == :constant
         cache.dtprev = one(cache.dtprev)
         uzero = zero(eltype(u))
-        for i in 1:s
+        for i in 1 : num_stages
             @..  z[i] = w[i] = uzero
         end
-        for i in 1:(s-1)
+        for i in 1 : (num_stages-1)
             @.. cache.cont[i] = uzero
         end
     else
-        c' = Vector{eltype(u)}(undef, s) #time stepping
-        c'[s] = dt / cache.dtprev
-        for i in 1 : s-1
-            c'[i] = c[i] * c'[s]
-        end
-        for i in 1 : s # collocation polynomial
-            @.. z[i] = cont[s-1] * (c'[i] - c[1] + 1) + cont[s-1]
-            j = s - 2
+        c' = Vector{eltype(u)}(undef, num_stages) #time stepping
+        c'[num_stages] = dt / cache.dtprev
+        for i in 1 : (num_stages-1)
+            c'[i] = c[i] * c'[num_stages]
+        end
+        for i in 1 : num_stages # collocation polynomial
+            @.. z[i] = cont[num_stages - 1] * (c'[i] - c[1] + 1) + cont[num_stages - 1]
+            j = num_stages - 2
             while j > 0
-                @.. z[i] = z[i] * (c'[i] - c[s-j] + 1) + cont[j]
+                @.. z[i] = z[i] * (c'[i] - c[num_stages - j] + 1) + cont[j]
             end
             @.. z[i] = z[i] * c'[i]
         end
@@ -1601,22 +1601,22 @@ end
 
         # evaluate function
         k[1] = fsallast
-        for i in 1 : s
+        for i in 1 : num_stages
             @.. tmp = uprev + z[i]
             f(k[i], tmp, p, t + c[i] * dt)
         end
-        integrator.stats.nf += s
+        integrator.stats.nf += num_stages
 
         @.. fw = TI * k
         if mass_matrix === I
             Mw = w
         elseif mass_matrix isa UniformScaling
-            for i in 1 : s
+            for i in 1 : num_stages
                 mul!(z[i], mass_matrix.λ, w[i])
             end
             Mw = z
         else
-            for i in 1 : s
+            for i in 1 : num_stages
                 mul!(z[i], mass_matrix.λ, w[i])
             end
             Mw = z
@@ -1626,7 +1626,7 @@ end
         needfactor = iter == 1 && new_W
 
         linsolve1 = cache.linsolve1
-        linres = Vector{BigFloat}(undef, s)
+        linres = Vector{BigFloat}(undef, num_stages)
         if needfactor
             linres[1] = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff),
                 linu = _vec(dw1))
@@ -1637,9 +1637,9 @@ end
 
         cache.linsolve1 = linres1.cache
 
-        for i in 1 : (s-1)/2
+        for i in 1 : (num_stages - 1)/2
             @.. broadcast=false cubuff[i]=complex(
-            fw2 - α[i]dt * Mw[2*i] + β[i]dt * Mw[2*i + 1], fw3 - β[i]dt * Mw[2*i] - α[i]dt * Mw[2*i + 1])
+            fw2 - αdt[i] * Mw[2 * i] + βdt[i] * Mw[2 * i + 1], fw3 - βdt[i] * Mw[2 * i] - αdt[i] * Mw[2 * i + 1])
             linsolve2[i] = cache.linsolve2[i]
             if needfactor
                 linres[i + 1] = dolinsolve(integrator, linsolve2[i]; A = W2[i], b = _vec(cubuff[i]),
@@ -1651,10 +1651,10 @@ end
             cache.linsolve2[i] = linres[i + 1].cache
         end
 
-        integrator.stats.nsolve += (s+1) / 2
+        integrator.stats.nsolve += (num_stages + 1) / 2
         dw[1] = dw1
         i = 2
-        while i <= s
+        while i <= num_stages
             dw[i] = z[i]
             dw[i + 1] = z[i + 1]
             @.. dw[i] = real(dw2[i - 1])
@@ -1664,14 +1664,14 @@ end
 
         # compute norm of residuals
         iter > 1 && (ndwprev = ndw)
-        ndws = Vector{BigFloat}(undef, s)
-        for i in 1:s
+        ndws = Vector{BigFloat}(undef, num_stages)
+        for i in 1:num_stages
             calculate_residuals!(atmp, dw[i], uprev, u, atol, rtol, internalnorm, t)
             ndws[i] = internalnorm(atmp, t)
         end
 
         ndw = 0
-        for i in 1 : s
+        for i in 1 : num_stages
             ndw += ndws[i]
         end
 
@@ -1714,7 +1714,7 @@ end
     @.. broadcast=false u=uprev + z[s]
 
     step_limiter!(u, integrator, p, t + dt)
-
+    #=
     if adaptive
         utilde = w2
         edt = e./dt
@@ -1752,13 +1752,13 @@ end
             integrator.EEst = internalnorm(atmp, t)
         end
     end
-
+    =#
     if integrator.EEst <= oneunit(integrator.EEst)
         cache.dtprev = dt
         if alg.extrapolant != :constant
-            derivatives = Matrix{eltype(u)}(undef, s-1, s-1)
-            for i in 1 : (s - 1)
-                for j in i : (s-1)
+            derivatives = Matrix{eltype(u)}(undef, num_stages - 1, num_stages - 1)
+            for i in 1 : (num_stages - 1)
+                for j in i : (num_stages - 1)
                     if i == 1
                         @.. derivatives[i, j] = (z[i] - z[i + 1]) / (c[i] - c[i + 1]) #first derivatives
                     else
@@ -1766,8 +1766,8 @@ end
                     end
                 end
             end
-            for i in 1 : (s-1)
-                cache.cont[i] = derivatives[i, s - 1]
+            for i in 1 : (num_stages - 1)
+                cache.cont[i] = derivatives[i, num_stages - 1]
             end
         end
     end
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index 478bdd4240..3e36fd9363 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -259,7 +259,7 @@ function RadauIIA9Tableau(T, T2)
         e1, e2, e3, e4, e5)
 end
 
-struct adaptiveRadauTableau(T, T2)
+struct adaptiveRadauTableau{T, T2, Int}
     T:: AbstractMatrix{T}
     TI::AbstractMatrix{T}
     γ::T
@@ -267,80 +267,81 @@ struct adaptiveRadauTableau(T, T2)
     β::AbstractVector{T}
     c::AbstractVector{T}
     e::AbstractVector{T}
+    S::Int
 end
 
 using Polynomials, GenericLinearAlgebra, LinearAlgebra, LinearSolve, GenericSchur, BSeries
 
-function adaptiveRadauTableau(T, T2, s::Int64)
-    tmp = Vector{BigFloat}(undef, s-1)
-    for i in 1:(s-1)
+function adaptiveRadauTableau(T, T2, num_stages::Int)
+    tmp = Vector{BigFloat}(undef, num_stages - 1)
+    for i in 1:(num_stages - 1)
         tmp[i] = 0
     end
-    tmp2 = Vector{BigFloat}(undef, s+1)
-    for i in 1:(s+1)
-       tmp2[i]=(-1)^(s+1-i) * binomial(s,s+1-i)
+    tmp2 = Vector{BigFloat}(undef, num_stages + 1)
+    for i in 1:(num_stages + 1)
+       tmp2[i]=(-1)^(num_stages + 1 - i) * binomial(num_stages , num_stages + 1 - i)
     end
     p = Polynomial{BigFloat}([tmp; tmp2])
-    for i in 1:s-1
+    for i in 1:(num_stages - 1)
         p = derivative(p)
     end
     c = roots(p)
-    c[s] = 1
-    c_powers = Matrix{BigFloat}(undef, s, s)
-    for i in 1:s
-        for j in 1:s
-            c_powers[i,j] = c[i]^(j-1)
+    c[num_stages] = 1
+    c_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
+    for i in 1 : num_stages
+        for j in 1 : num_stages
+            c_powers[i,j] = c[i]^(j - 1)
         end
     end
     inverse_c_powers = inv(c_powers)
-    c_q = Matrix{BigFloat}(undef, s, s)
-    for i in 1:s
-        for j in 1:s
+    c_q = Matrix{BigFloat}(undef, num_stages, num_stages)
+    for i in 1 : num_stages
+        for j in 1 : num_stages
             c_q[i,j] = c[i]^(j) / j
         end
     end
     a = c_q * inverse_c_powers
     a_inverse = inv(a)
-    b = Vector{Bigfloat}(undef, s)
-    for i in 1 : s
-        b[i] = a[s, i]
+    b = Vector{BigFloat}(undef, num_stages)
+    for i in 1 : num_stages
+        b[i] = a[num_stages, i]
     end
     vals = eigvals(a_inverse)
-    γ = real(b[s])
-    α = Vector{BigFloat}(undef, floor(Int, s/2))
-    β = Vector{BigFloat}(undef, floor(Int, s/2))
+    γ = real(b[num_stages])
+    α = Vector{BigFloat}(undef, floor(Int, num_stages/2))
+    β = Vector{BigFloat}(undef, floor(Int, num_stages/2))
     index = 1
     i = 1
-    while i <= (s-1)
+    while i <= (num_stages - 1)
         α[index] = real(vals[i])
         β[index] = imag(vals[i + 1])
         index = index + 1
         i = i + 2
     end
     eigvec = eigvecs(a)
-    vecs = Vector{Vector{BigFloat}}(undef, s)
+    vecs = Vector{Vector{BigFloat}}(undef, num_stages)
     i = 1
     index = 2
-    while i < s 
-        vecs[index] = real(eigvec[:, i] ./ eigvec[s,i])
-        vecs[index + 1] = -imag(eigvec[:, i] ./ eigvec[s,i])
+    while i < num_stages 
+        vecs[index] = real(eigvec[:, i] ./ eigvec[num_stages, i])
+        vecs[index + 1] = -imag(eigvec[:, i] ./ eigvec[num_stages, i])
         index += 2
         i += 2
     end
-    vecs[1] = real(eigvec[:, s])
+    vecs[1] = real(eigvec[:, num_stages])
     tmp3 = vcat(vecs)
-    T = Matrix{BigFloat}(undef, s, s)
-    for j in 1 : s
-        for i in 1 : s
+    T = Matrix{BigFloat}(undef, num_stages, num_stages)
+    for j in 1 : num_stages
+        for i in 1 : num_stages
             T[i, j] = tmp3[j][i]
         end
     end
     TI = inv(T)
-    b_hat = Vector{BigFloat}(undef, s)
-    embedded = bseries(a, b_hat, c, s - 2)
+    #b_hat = Vector{BigFloat}(undef, num_stages)
+    #embedded = bseries(a, b_hat, c, num_stages - 2)
 
     #e = b_hat - b
-    #adaptiveRadautableau(T, TI, γ, α, β, c, e)
+    #adaptiveRadautableau(T, TI, γ, α, β, c, e, s)
 end
 
 adaptiveRadauTableau(0, 0, 3)
diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index 5cc4bfc597..b68f3cc72e 100644
--- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -8,6 +8,8 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear]
     @test sim21.𝒪est[:final]≈5 atol=testTol
 end
 
+sol = solve(prob_ode_linear, AdaptiveRadau())
+
 sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9())
 @test sim21.𝒪est[:final]≈8 atol=testTol