diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 78cb0d9372..296ae17352 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -90,8 +90,8 @@ jobs:
         env:
           GROUP: ${{ matrix.group }}
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info
-          fail_ci_if_error: true
+          fail_ci_if_error: false
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
index 92c0519266..a8d5a89a9a 100644
--- a/.github/workflows/CompatHelper.yml
+++ b/.github/workflows/CompatHelper.yml
@@ -23,4 +23,4 @@ jobs:
       - name: CompatHelper.main()
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-        run: julia -e 'using CompatHelper; CompatHelper.main(;subdirs=["", "docs", "test/downstream"])'
+        run: julia -e 'using CompatHelper; CompatHelper.main(;subdirs=["", "docs", "test/downstream", "lib"])'
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index 70c5a8e0ac..accac77277 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -32,7 +32,7 @@ jobs:
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
         run: julia --project=docs/ --code-coverage=user docs/make.jl
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info
diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml
index 8200ba42ba..dd57d15a64 100644
--- a/.github/workflows/Downstream.yml
+++ b/.github/workflows/Downstream.yml
@@ -30,7 +30,10 @@ jobs:
           - {user: SciML, repo: SciMLSensitivity.jl, group: Core3}
           - {user: SciML, repo: SciMLSensitivity.jl, group: Core4}
           - {user: SciML, repo: SciMLSensitivity.jl, group: Core5}
-          - {user: SciML, repo: ModelingToolkit.jl, group: All}
+          - {user: SciML, repo: ModelingToolkit.jl, group: InterfaceI}
+          - {user: SciML, repo: ModelingToolkit.jl, group: InterfaceII}
+          - {user: SciML, repo: ModelingToolkit.jl, group: Initialization}
+          - {user: SciML, repo: ModelingToolkit.jl, group: SymbolicIndexingInterface}
           - {user: SciML, repo: DiffEqDevTools.jl, group: Core}
           - {user: nathanaelbosch, repo: ProbNumDiffEq.jl, group: Downstream}
           - {user: SKopecz, repo: PositiveIntegrators.jl, group: Downstream}
@@ -73,7 +76,7 @@ jobs:
             exit(0)  # Exit immediately, as a success
           end
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info
diff --git a/Project.toml b/Project.toml
index da47ef934b..bd6935f76f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEq"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <mayingbo5@gmail.com>"]
-version = "6.89.0"
+version = "6.90.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml
index 3ebfcd37f1..f5f381c9a0 100644
--- a/lib/OrdinaryDiffEqCore/Project.toml
+++ b/lib/OrdinaryDiffEqCore/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqCore"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "1.10.0"
+version = "1.13.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -49,7 +49,7 @@ Accessors = "0.1.36"
 Adapt = "3.0, 4"
 ArrayInterface = "7"
 DataStructures = "0.18"
-DiffEqBase = "6.157"
+DiffEqBase = "6.160"
 DiffEqDevTools = "2.44.4"
 DocStringExtensions = "0.9"
 EnumX = "1"
@@ -70,7 +70,7 @@ Random = "<0.0.1, 1"
 RecursiveArrayTools = "2.36, 3"
 Reexport = "1.0"
 SafeTestsets = "0.1.0"
-SciMLBase = "2.57.2"
+SciMLBase = "2.62"
 SciMLOperators = "0.3"
 SciMLStructures = "1"
 SimpleUnPack = "1"
diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
index a7cbd95167..fbfcd66dd3 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -60,9 +60,9 @@ using DiffEqBase: check_error!, @def, _vec, _reshape
 
 using FastBroadcast: @.., True, False
 
-using SciMLBase: NoInit, CheckInit, _unwrap_val
+using SciMLBase: NoInit, CheckInit, OverrideInit, AbstractDEProblem, _unwrap_val
 
-import SciMLBase: alg_order
+import SciMLBase: AbstractNonlinearProblem, alg_order
 
 import DiffEqBase: calculate_residuals,
                    calculate_residuals!, unwrap_cache,
@@ -76,7 +76,8 @@ import Accessors: @reset
 
 using SciMLStructures: canonicalize, Tunable, isscimlstructure
 
-using SymbolicIndexingInterface: parameter_values, is_variable, variable_index, symbolic_type, NotSymbolic
+using SymbolicIndexingInterface: state_values, parameter_values, is_variable, variable_index,
+                                 symbolic_type, NotSymbolic
 
 const CompiledFloats = Union{Float32, Float64}
 import Preferences
diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl
index 5d144fcfae..d5f89349d6 100644
--- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl
@@ -17,6 +17,9 @@ end
 
 SciMLBase.forwarddiffs_model_time(alg::RosenbrockAlgorithm) = true
 
+SciMLBase.allows_late_binding_tstops(::OrdinaryDiffEqAlgorithm) = true
+SciMLBase.allows_late_binding_tstops(::DAEAlgorithm) = true
+
 # isadaptive is defined below.
 
 ## OrdinaryDiffEq Internal Traits
diff --git a/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl b/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl
index f836d61b9b..72f9a6cdfd 100644
--- a/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl
+++ b/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl
@@ -73,7 +73,8 @@ function alg_cache(alg::CompositeAlgorithm{CS, Tuple{A1, A2, A3, A4, A5, A6}}, u
     args = (u, rate_prototype, uEltypeNoUnits,
         uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt,
         reltol, p, calck, Val(V))
-    argT = map(typeof, args)
+    # Core.Typeof to turn uEltypeNoUnits into Type{uEltypeNoUnits} rather than DataType
+    argT = map(Core.Typeof, args)
     T1 = Base.promote_op(alg_cache, A1, argT...)
     T2 = Base.promote_op(alg_cache, A2, argT...)
     T3 = Base.promote_op(alg_cache, A3, argT...)
diff --git a/lib/OrdinaryDiffEqCore/src/dense/generic_dense.jl b/lib/OrdinaryDiffEqCore/src/dense/generic_dense.jl
index f8dfa068e2..9978478560 100644
--- a/lib/OrdinaryDiffEqCore/src/dense/generic_dense.jl
+++ b/lib/OrdinaryDiffEqCore/src/dense/generic_dense.jl
@@ -168,10 +168,12 @@ function default_ode_interpolant(
         return ode_interpolant(Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache5, idxs,
             deriv, integrator.differential_vars)
-    else # alg_choice == 6
+    elseif alg_choice == 6
         return ode_interpolant(Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache6, idxs,
             deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -227,6 +229,8 @@ end
             ode_interpolant!(val, Θ, integrator.dt, integrator.uprev, integrator.u,
                 integrator.k, integrator.cache.cache6,
                 idxs, deriv, integrator.differential_vars)
+        else
+            error("DefaultCache invalid alg_choice. File an issue.")
         end
     else
         ode_interpolant!(val, Θ, integrator.dt, integrator.uprev, integrator.u,
@@ -256,10 +260,12 @@ function default_ode_interpolant!(
         return ode_interpolant!(val, Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache5, idxs,
             deriv, integrator.differential_vars)
-    else # alg_choice == 6
+    elseif alg_choice == 6
         return ode_interpolant!(val, Θ, integrator.dt, integrator.uprev,
             integrator.u, integrator.k, cache.cache6, idxs,
             deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -380,6 +386,8 @@ function default_ode_extrapolant!(
         ode_interpolant!(val, Θ, integrator.t - integrator.tprev,
             integrator.uprev2, integrator.uprev,
             integrator.k, cache.cache6, idxs, deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -444,6 +452,8 @@ function default_ode_extrapolant(
         ode_interpolant(Θ, integrator.t - integrator.tprev,
             integrator.uprev2, integrator.uprev,
             integrator.k, cache.cache6, idxs, deriv, integrator.differential_vars)
+    else
+        error("DefaultCache invalid alg_choice. File an issue.")
     end
 end
 
@@ -810,6 +820,8 @@ function ode_interpolation(tval::Number, id::I, idxs, deriv::D, p,
                     cache.cache6) # update the kcurrent
                 val = ode_interpolant(Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊],
                     cache.cache6, idxs, deriv, differential_vars)
+            else
+                error("DefaultCache invalid alg_choice. File an issue.")
             end
         else
             _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p,
@@ -892,6 +904,8 @@ function ode_interpolation!(out, tval::Number, id::I, idxs, deriv::D, p,
                     cache.cache6) # update the kcurrent
                 ode_interpolant!(out, Θ, dt, timeseries[i₋], timeseries[i₊], ks[i₊],
                     cache.cache6, idxs, deriv, differential_vars)
+            else
+                error("DefaultCache invalid alg_choice. File an issue.")
             end
         else
             _ode_addsteps!(ks[i₊], ts[i₋], timeseries[i₋], timeseries[i₊], dt, f, p,
diff --git a/lib/OrdinaryDiffEqCore/src/initialize_dae.jl b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl
index ccb2db0a7f..4e156175c0 100644
--- a/lib/OrdinaryDiffEqCore/src/initialize_dae.jl
+++ b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl
@@ -20,16 +20,6 @@ function BrownFullBasicInit(; abstol = 1e-10, nlsolve = nothing)
 end
 BrownFullBasicInit(abstol) = BrownFullBasicInit(; abstol = abstol, nlsolve = nothing)
 
-struct OverrideInit{T, F} <: DiffEqBase.DAEInitializationAlgorithm
-    abstol::T
-    nlsolve::F
-end
-
-function OverrideInit(; abstol = 1e-10, nlsolve = nothing)
-    OverrideInit(abstol, nlsolve)
-end
-OverrideInit(abstol) = OverrideInit(; abstol = abstol, nlsolve = nothing)
-
 ## Notes
 
 #=
@@ -116,7 +106,7 @@ default_nlsolve(alg, isinplace, u, initprob, autodiff = false) = alg
 
 ## If the initialization is trivial just use nothing alg
 function default_nlsolve(
-        ::Nothing, isinplace::Val{true}, u::Nothing, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace::Val{true}, u::Nothing, ::AbstractNonlinearProblem, autodiff = false)
     nothing
 end
 
@@ -126,7 +116,7 @@ function default_nlsolve(
 end
 
 function default_nlsolve(
-        ::Nothing, isinplace::Val{false}, u::Nothing, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace::Val{false}, u::Nothing, ::AbstractNonlinearProblem, autodiff = false)
     nothing
 end
 
@@ -137,7 +127,7 @@ function default_nlsolve(
 end
 
 function OrdinaryDiffEqCore.default_nlsolve(
-        ::Nothing, isinplace, u, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace, u, ::AbstractNonlinearProblem, autodiff = false)
     error("This ODE requires a DAE initialization and thus a nonlinear solve but no nonlinear solve has been loaded. To solve this problem, do `using OrdinaryDiffEqNonlinearSolve` or pass a custom `nlsolve` choice into the `initializealg`.")
 end
 
@@ -148,24 +138,21 @@ end
 
 ## NoInit
 
-function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
+function _initialize_dae!(integrator, prob::AbstractDEProblem,
         alg::NoInit, x::Union{Val{true}, Val{false}})
 end
 
 ## OverrideInit
 
-function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
+function _initialize_dae!(integrator, prob::AbstractDEProblem,
         alg::OverrideInit, isinplace::Union{Val{true}, Val{false}})
-    initializeprob = prob.f.initializeprob
-
-    if SciMLBase.has_update_initializeprob!(prob.f)
-        prob.f.update_initializeprob!(initializeprob, prob)
-    end
+    initializeprob = prob.f.initialization_data.initializeprob
 
     # If it doesn't have autodiff, assume it comes from symbolic system like ModelingToolkit
     # Since then it's the case of not a DAE but has initializeprob
     # In which case, it should be differentiable
-    isAD = if initializeprob.u0 === nothing
+    iu0 = state_values(initializeprob)
+    isAD = if iu0 === nothing
         AutoForwardDiff
     elseif has_autodiff(integrator.alg)
         alg_autodiff(integrator.alg) isa AutoForwardDiff
@@ -173,115 +160,30 @@ function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
         true
     end
 
-    alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD)
-    nlsol = solve(initializeprob, alg)
+    nlsolve_alg = default_nlsolve(alg.nlsolve, isinplace, iu0, initializeprob, isAD)
+
+    u0, p, success = SciMLBase.get_initial_values(prob, integrator, prob.f, alg, isinplace; nlsolve_alg, abstol = integrator.opts.abstol, reltol = integrator.opts.reltol)
+
     if isinplace === Val{true}()
-        integrator.u .= prob.f.initializeprobmap(nlsol)
+        integrator.u .= u0
     elseif isinplace === Val{false}()
-        integrator.u = prob.f.initializeprobmap(nlsol)
+        integrator.u = u0
     else
         error("Unreachable reached. Report this error.")
     end
-    if SciMLBase.has_initializeprobpmap(prob.f)
-        integrator.p = prob.f.initializeprobpmap(prob, nlsol)
-        sol = integrator.sol
-        @reset sol.prob.p = integrator.p
-        integrator.sol = sol
-    end
+    integrator.p = p
+    sol = integrator.sol
+    @reset sol.prob.p = integrator.p
+    integrator.sol = sol
 
-    if nlsol.retcode != ReturnCode.Success
+    if !success
         integrator.sol = SciMLBase.solution_new_retcode(integrator.sol,
             ReturnCode.InitialFailure)
     end
 end
 
 ## CheckInit
-struct CheckInitFailureError <: Exception
-    normresid::Any
-    abstol::Any
-end
-
-function Base.showerror(io::IO, e::CheckInitFailureError)
-    print(io,
-        "DAE initialization failed: your u0 did not satisfy the initialization requirements, 
-        normresid = $(e.normresid) > abstol = $(e.abstol). If you wish for the system to 
-        automatically change the algebraic variables to satisfy the algebraic constraints, 
-        please pass `initializealg = BrownBasicInit()` to solve (this option will require 
-        `using OrdinaryDiffEqNonlinearSolve`). If you wish to perform an initialization on the
-        complete u0, please pass initializealg = ShampineCollocationInit() to solve. Note that 
-        initialization can be a very difficult process for DAEs and in many cases can be 
-        numerically intractable without symbolic manipulation of the system. For an automated 
-        system that will generate numerically stable initializations, see ModelingToolkit.jl 
-        structural simplification for more details."
-    )
-end
-
-function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit,
-        isinplace::Val{true})
-    @unpack p, t, f = integrator
-    M = integrator.f.mass_matrix
-    tmp = first(get_tmp_cache(integrator))
-    u0 = integrator.u
-
-    algebraic_vars = [all(iszero, x) for x in eachcol(M)]
-    algebraic_eqs = [all(iszero, x) for x in eachrow(M)]
-    (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
-    update_coefficients!(M, u0, p, t)
-    f(tmp, u0, p, t)
-    tmp .= ArrayInterface.restructure(tmp, algebraic_eqs .* _vec(tmp))
-
-    normresid = integrator.opts.internalnorm(tmp, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
-end
-
-function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit,
-        isinplace::Val{false})
-    @unpack p, t, f = integrator
-    u0 = integrator.u
-    M = integrator.f.mass_matrix
-
-    algebraic_vars = [all(iszero, x) for x in eachcol(M)]
-    algebraic_eqs = [all(iszero, x) for x in eachrow(M)]
-    (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
-    update_coefficients!(M, u0, p, t)
-    du = f(u0, p, t)
-    resid = _vec(du)[algebraic_eqs]
-
-    normresid = integrator.opts.internalnorm(resid, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
-end
-
-function _initialize_dae!(integrator, prob::DAEProblem,
-        alg::CheckInit, isinplace::Val{true})
-    @unpack p, t, f = integrator
-    u0 = integrator.u
-    resid = get_tmp_cache(integrator)[2]
-
-    f(resid, integrator.du, u0, p, t)
-    normresid = integrator.opts.internalnorm(resid, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
-end
-
-function _initialize_dae!(integrator, prob::DAEProblem,
-        alg::CheckInit, isinplace::Val{false})
-    @unpack p, t, f = integrator
-    u0 = integrator.u
-
-    nlequation_oop = u -> begin
-        f((u - u0) / dt, u, p, t)
-    end
-
-    nlequation = (u, _) -> nlequation_oop(u)
-
-    resid = f(integrator.du, u0, p, t)
-    normresid = integrator.opts.internalnorm(resid, t)
-    if normresid > integrator.opts.abstol
-        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
-    end
+function _initialize_dae!(integrator, prob::AbstractDEProblem, alg::CheckInit,
+        isinplace::Union{Val{true}, Val{false}})
+    SciMLBase.get_initial_values(prob, integrator, prob.f, alg, isinplace; abstol = integrator.opts.abstol)
 end
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
index 9a0bb53267..a38ce1566a 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
@@ -400,3 +400,56 @@ function post_newton_controller!(integrator, alg)
     integrator.dt = integrator.dt / integrator.opts.failfactor
     nothing
 end
+
+@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
+    @unpack qmin, qmax, gamma = integrator.opts
+    EEst = DiffEqBase.value(integrator.EEst)
+    if iszero(EEst)
+        q = inv(qmax)
+    else
+        if fac_default_gamma(alg)
+            fac = gamma
+        else
+            if isfirk(alg)
+                @unpack iter = integrator.cache
+                @unpack maxiters = alg
+            else
+                @unpack iter, maxiters = integrator.cache.nlsolver
+            end
+            fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
+        end
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qtmp = FastPower.fastpower(EEst, expo) / fac
+        @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
+        integrator.qold = q
+    end
+    q
+end
+
+function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
+    @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
+ 
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if integrator.success_iter > 0
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qgus = (integrator.dtacc / integrator.dt) *
+               FastPower.fastpower((EEst^2) / integrator.erracc, expo)
+        qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
+        qacc = max(q, qgus)
+    else
+        qacc = q
+    end
+    if qsteady_min <= qacc <= qsteady_max
+        qacc = one(qacc)
+    end
+    integrator.dtacc = integrator.dt
+    integrator.erracc = max(1e-2, EEst)
+    
+    return integrator.dt / qacc
+end
+
+function step_reject_controller!(integrator, controller::PredictiveController, alg)
+    @unpack dt, success_iter, qold = integrator
+    integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
+end
\ No newline at end of file
diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl
index ef77e4c6b4..42fd07a5bf 100644
--- a/lib/OrdinaryDiffEqCore/src/solve.jl
+++ b/lib/OrdinaryDiffEqCore/src/solve.jl
@@ -242,6 +242,12 @@ function DiffEqBase.__init(
         resType = typeof(res_prototype)
     end
 
+    if tstops isa AbstractArray || tstops isa Tuple || tstops isa Number
+        _tstops = nothing
+    else
+        _tstops = tstops
+        tstops = ()
+    end
     tstops_internal = initialize_tstops(tType, tstops, d_discontinuities, tspan)
     saveat_internal = initialize_saveat(tType, saveat, tspan)
     d_discontinuities_internal = initialize_d_discontinuities(tType, d_discontinuities,
@@ -264,29 +270,7 @@ function DiffEqBase.__init(
     end
 
     ### Algorithm-specific defaults ###
-    if save_idxs === nothing
-        saved_subsystem = nothing
-    else
-        if !(save_idxs isa AbstractArray) || symbolic_type(save_idxs) != NotSymbolic()
-            _save_idxs = [save_idxs]
-        else
-            _save_idxs = save_idxs
-        end
-        saved_subsystem = SciMLBase.SavedSubsystem(prob, parameter_values(prob), _save_idxs)
-        if saved_subsystem !== nothing
-            _save_idxs = SciMLBase.get_saved_state_idxs(saved_subsystem)
-            if isempty(_save_idxs)
-                # no states to save
-                save_idxs = Int[]
-            elseif !(save_idxs isa AbstractArray) || symbolic_type(save_idxs) != NotSymbolic()
-                # only a single state to save, and save it as a scalar timeseries instead of
-                # single-element array
-                save_idxs = only(_save_idxs)
-            else
-                save_idxs = _save_idxs
-            end
-        end
-    end
+    save_idxs, saved_subsystem = SciMLBase.get_save_idxs_and_saved_subsystem(prob, save_idxs)
 
     if save_idxs === nothing
         ksEltype = Vector{rateType}
@@ -564,6 +548,13 @@ function DiffEqBase.__init(
         end
     end
 
+    if _tstops !== nothing
+        tstops = _tstops(parameter_values(integrator), prob.tspan)
+        for tstop in tstops
+            add_tstop!(integrator, tstop)
+        end
+    end
+
     handle_dt!(integrator)
     integrator
 end
diff --git a/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl b/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl
index 8676b4a57a..c434bc80ab 100644
--- a/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl
+++ b/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl
@@ -40,6 +40,8 @@ end
 prob_rober = ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1e3), (0.04, 3e7, 1e4))
 sol = solve(prob_rober)
 rosensol = solve(prob_rober, AutoTsit5(Rosenbrock23(autodiff = false)))
+#test that cache is type stable
+@test typeof(sol.interp.cache.cache3) == typeof(rosensol.interp.cache.caches[2])
 # test that default has the same performance as AutoTsit5(Rosenbrock23()) (which we expect it to use for this).
 @test sol.stats.naccept == rosensol.stats.naccept
 @test sol.stats.nf == rosensol.stats.nf
@@ -50,6 +52,8 @@ rosensol = solve(prob_rober, AutoTsit5(Rosenbrock23(autodiff = false)))
 sol = solve(prob_rober, reltol = 1e-7, abstol = 1e-7)
 rosensol = solve(
     prob_rober, AutoVern7(Rodas5P(autodiff = false)), reltol = 1e-7, abstol = 1e-7)
+#test that cache is type stable
+@test typeof(sol.interp.cache.cache4) == typeof(rosensol.interp.cache.caches[2])
 # test that default has the same performance as AutoTsit5(Rosenbrock23()) (which we expect it to use for this).
 @test sol.stats.naccept == rosensol.stats.naccept
 @test sol.stats.nf == rosensol.stats.nf
diff --git a/lib/OrdinaryDiffEqDifferentiation/Project.toml b/lib/OrdinaryDiffEqDifferentiation/Project.toml
index f1a8b2371c..11d8b4cfb0 100644
--- a/lib/OrdinaryDiffEqDifferentiation/Project.toml
+++ b/lib/OrdinaryDiffEqDifferentiation/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqDifferentiation"
 uuid = "4302a76b-040a-498a-8c04-15b101fed76b"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <mayingbo5@gmail.com>"]
-version = "1.1.0"
+version = "1.2.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml
index a0cc91007c..62fc014d3c 100644
--- a/lib/OrdinaryDiffEqFIRK/Project.toml
+++ b/lib/OrdinaryDiffEqFIRK/Project.toml
@@ -1,34 +1,30 @@
 name = "OrdinaryDiffEqFIRK"
 uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "1.2.0"
+version = "1.5.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
-GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
-GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b"
 OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
-Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
-RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 DiffEqBase = "6.152.2"
 DiffEqDevTools = "2.44.4"
 FastBroadcast = "0.3.5"
+FastGaussQuadrature = "1.0.2"
 FastPower = "1"
-GenericLinearAlgebra = "0.3.13"
-GenericSchur = "0.5.4"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"
 MuladdMacro = "0.2.4"
@@ -36,23 +32,22 @@ ODEProblemLibrary = "0.1.8"
 OrdinaryDiffEqCore = "1.1"
 OrdinaryDiffEqDifferentiation = "<0.0.1, 1"
 OrdinaryDiffEqNonlinearSolve = "<0.0.1, 1"
-Polynomials = "4.0.11"
 Random = "<0.0.1, 1"
 RecursiveArrayTools = "3.27.0"
 Reexport = "1.2.2"
-RootedTrees = "2.23.1"
 SafeTestsets = "0.1.0"
+SciMLBase = "2.60.0"
 SciMLOperators = "0.3.9"
-Symbolics = "6.15.3"
 Test = "<0.0.1, 1"
 julia = "1.10"
 
 [extras]
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test", "ODEProblemLibrary"]
+test = ["DiffEqDevTools", "GenericSchur", "Random", "SafeTestsets", "Test", "ODEProblemLibrary"]
diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
index 753f094704..5817abd9b7 100644
--- a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
@@ -18,7 +18,6 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            get_current_adaptive_order, get_fsalfirstlast,
                            isfirk, generic_solver_docstring
 using MuladdMacro, DiffEqBase, RecursiveArrayTools
-using Polynomials, GenericLinearAlgebra, GenericSchur
 using SciMLOperators: AbstractSciMLOperator
 using LinearAlgebra: I, UniformScaling, mul!, lu
 import LinearSolve
diff --git a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
index 428393e0bf..ef91903e63 100644
--- a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
@@ -1,9 +1,9 @@
-qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}) = 8
+qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA9, AdaptiveRadau}) = 8
 
 alg_order(alg::RadauIIA3) = 3
 alg_order(alg::RadauIIA5) = 5
 alg_order(alg::RadauIIA9) = 9
-alg_order(alg::AdaptiveRadau) = 5
+alg_order(alg::AdaptiveRadau) = 5 #dummy value
 
 isfirk(alg::RadauIIA3) = true
 isfirk(alg::RadauIIA5) = true
@@ -13,3 +13,6 @@ isfirk(alg::AdaptiveRadau) = true
 alg_adaptive_order(alg::RadauIIA3) = 1
 alg_adaptive_order(alg::RadauIIA5) = 3
 alg_adaptive_order(alg::RadauIIA9) = 5
+
+get_current_alg_order(alg::AdaptiveRadau, cache) = cache.num_stages * 2 - 1
+get_current_adaptive_order(alg::AdaptiveRadau, cache) = cache.num_stages
diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
index 25165422c8..e98389ec7c 100644
--- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -163,12 +163,13 @@ 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
+ min_order::Int
+ max_order::Int
 end
 
 function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      standardtag = Val{true}(), concrete_jac = nothing,
-     diff_type = Val{:forward}, num_stages = 3, 
+     diff_type = Val{:forward}, min_order = 5, max_order = 13, 
      linsolve = nothing, precs = DEFAULT_PRECS,
      extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
      new_W_γdt_cutoff = 1 // 5,
@@ -186,6 +187,6 @@ function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      fast_convergence_cutoff,
      new_W_γdt_cutoff,
      controller,
-     step_limiter!, num_stages)
+     step_limiter!, min_order, max_order)
 end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
index 2887bacd4f..3849d8114c 100644
--- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -1,37 +1,14 @@
-@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
-    @unpack qmin, qmax, gamma = integrator.opts
-    EEst = DiffEqBase.value(integrator.EEst)
-
-    if iszero(EEst)
-        q = inv(qmax)
-    else
-        if fac_default_gamma(alg)
-            fac = gamma
-        else
-            if isfirk(alg)
-                @unpack iter = integrator.cache
-                @unpack maxiters = alg
-            else
-                @unpack iter, maxiters = integrator.cache.nlsolver
-            end
-            fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
-        end
-        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
-        qtmp = FastPower.fastpower(EEst, expo) / fac
-        @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
-        integrator.qold = q
-    end
-    q
-end
-
-function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
+function step_accept_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau, q)
     @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
+    @unpack cache = integrator
+    @unpack num_stages, step, iter, hist_iter = cache
+ 
     EEst = DiffEqBase.value(integrator.EEst)
 
     if integrator.success_iter > 0
         expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
         qgus = (integrator.dtacc / integrator.dt) *
-               FastPower.fastpower((EEst^2) / integrator.erracc, expo)
+               DiffEqBase.fastpow((EEst^2) / integrator.erracc, expo)
         qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
         qacc = max(q, qgus)
     else
@@ -42,10 +19,40 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
     end
     integrator.dtacc = integrator.dt
     integrator.erracc = max(1e-2, EEst)
+    cache.step = step + 1
+    hist_iter = hist_iter * 0.8 + iter * 0.2
+    cache.hist_iter = hist_iter
+    max_stages = (alg.max_order - 1) ÷ 4 * 2 + 1
+    min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1
+    if (step > 10)
+        if (hist_iter < 2.6 && num_stages <= max_stages)
+            cache.num_stages += 2
+            cache.step = 1
+            cache.hist_iter = iter
+        elseif ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages >= min_stages)
+            cache.num_stages -= 2 
+            cache.step = 1
+            cache.hist_iter = iter
+        end
+    end
     return integrator.dt / qacc
 end
 
-function step_reject_controller!(integrator, controller::PredictiveController, alg)
+function step_reject_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau)
     @unpack dt, success_iter, qold = integrator
+    @unpack cache = integrator 
+    @unpack num_stages, step, iter, hist_iter = cache
     integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
+    cache.step = step + 1
+    hist_iter = hist_iter * 0.8 + iter * 0.2
+    cache.hist_iter = hist_iter
+    min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1
+    if (step > 10)
+        if ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages >= min_stages)
+            cache.num_stages -= 2 
+            cache.step = 1
+            cache.hist_iter = iter
+        end
+    end
 end
+
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index af6da2a282..fb773cf4b0 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -362,6 +362,10 @@ mutable struct RadauIIA9Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     tmp4::uType
     tmp5::uType
     tmp6::uType
+    tmp7::uType
+    tmp8::uType
+    tmp9::uType
+    tmp10::uType
     atmp::uNoUnitsType
     jac_config::JC
     linsolve1::F1
@@ -440,6 +444,10 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits},
     tmp4 = zero(u)
     tmp5 = zero(u)
     tmp6 = zero(u)
+    tmp7 = zero(u)
+    tmp8 = zero(u)
+    tmp9 = zero(u)
+    tmp10 = zero(u)
     atmp = similar(u, uEltypeNoUnits)
     recursivefill!(atmp, false)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1)
@@ -469,7 +477,7 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits},
         du1, fsalfirst, k, k2, k3, k4, k5, fw1, fw2, fw3, fw4, fw5,
         J, W1, W2, W3,
         uf, tab, κ, one(uToltype), 10000,
-        tmp, tmp2, tmp3, tmp4, tmp5, tmp6, atmp, jac_config,
+        tmp, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, atmp, jac_config,
         linsolve1, linsolve2, linsolve3, rtol, atol, dt, dt,
         Convergence, alg.step_limiter!)
 end
@@ -477,7 +485,7 @@ end
 mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <:
     OrdinaryDiffEqConstantCache
     uf::F
-    tab::Tab
+    tabs::Vector{Tab}
     κ::Tol
     ηold::Tol
     iter::Int
@@ -486,6 +494,9 @@ mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <:
     W_γdt::Dt
     status::NLStatus
     J::JType
+    num_stages::Int
+    step::Int
+    hist_iter::Float64
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -494,30 +505,33 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     uf = UDerivativeWrapper(f, t, p)
     uToltype = constvalue(uBottomEltypeNoUnits)
-    num_stages = alg.num_stages
-
-    if (num_stages == 3)
-        tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 5)
-        tab = BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 7)
-        tab = BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif iseven(num_stages) || num_stages <3
-        error("num_stages must be odd and 3 or greater")
-    else
-        tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages)
+
+    max_order = alg.max_order
+    min_order = alg.min_order
+    max = (max_order - 1) ÷ 4 * 2 + 1
+    min = (min_order - 1) ÷ 4 * 2 + 1
+    if (alg.min_order < 5)
+        error("min_order choice $min_order below 5 is not compatible with the algorithm")
+    elseif (max < min)
+        error("max_order $max_order is below min_order $min_order")
     end
+    num_stages = min
 
-    cont = Vector{typeof(u)}(undef, num_stages)
-    for i in 1: num_stages
+    tabs = [RadauIIATableau5(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau9(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau13(uToltype, constvalue(tTypeNoUnits))]
+    i = 9
+    while i <= max
+        push!(tabs, RadauIIATableau(uToltype, constvalue(tTypeNoUnits), i))
+        i += 2
+    end
+    cont = Vector{typeof(u)}(undef, max)
+    for i in 1:max
         cont[i] = zero(u)
     end
 
     κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
     J = false .* _vec(rate_prototype) .* _vec(rate_prototype)'
-
-    AdaptiveRadauConstantCache(uf, tab, κ, one(uToltype), 10000, cont, dt, dt,
-        Convergence, J)
+    AdaptiveRadauConstantCache(uf, tabs, κ, one(uToltype), 10000, cont, dt, dt,
+        Convergence, J, num_stages, 1, 0.0)
 end
 
 mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, JType, W1Type, W2Type,
@@ -528,6 +542,8 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     z::Vector{uType}
     w::Vector{uType}
     c_prime::Vector{tType}
+    αdt::Vector{tType}
+    βdt::Vector{tType}
     dw1::uType
     ubuff::uType
     dw2::Vector{cuType}
@@ -544,7 +560,7 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     W1::W1Type #real
     W2::Vector{W2Type} #complex
     uf::UF
-    tab::Tab
+    tabs::Vector{Tab}
     κ::Tol
     ηold::Tol
     iter::Int
@@ -559,6 +575,9 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     W_γdt::Dt
     status::NLStatus
     step_limiter!::StepLimiter
+    num_stages::Int
+    step::Int
+    hist_iter::Float64
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -567,54 +586,60 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     uf = UJacobianWrapper(f, t, p)
     uToltype = constvalue(uBottomEltypeNoUnits)
-    num_stages = alg.num_stages
-
-    if (num_stages == 3)
-        tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 5)
-        tab = BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 7)
-        tab = BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif iseven(num_stages) || num_stages < 3
-        error("num_stages must be odd and 3 or greater")
-    else
-        tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages)
+
+    max_order = alg.max_order
+    min_order = alg.min_order
+    max = (max_order - 1) ÷ 4 * 2 + 1
+    min = (min_order - 1) ÷ 4 * 2 + 1
+    if (alg.min_order < 5)
+        error("min_order choice $min_order below 5 is not compatible with the algorithm")
+    elseif (max < min)
+        error("max_order $max_order is below min_order $min_order")
+    end
+    num_stages = min
+
+    tabs = [RadauIIATableau5(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau9(uToltype, constvalue(tTypeNoUnits)), RadauIIATableau13(uToltype, constvalue(tTypeNoUnits))]
+    i = 9
+    while i <= max
+        push!(tabs, RadauIIATableau(uToltype, constvalue(tTypeNoUnits), i))
+        i += 2
     end
 
     κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
 
-    z = Vector{typeof(u)}(undef, num_stages)
-    w = Vector{typeof(u)}(undef, num_stages)
-    for i in 1 : num_stages
-        z[i] = w[i] = zero(u)
+    z = Vector{typeof(u)}(undef, max)
+    w = Vector{typeof(u)}(undef, max)
+    for i in 1 : max
+        z[i] = zero(u)
+        w[i] = zero(u)
     end
 
-    c_prime = Vector{typeof(t)}(undef, num_stages) #time stepping
+    αdt = [zero(t) for i in 1:max]
+    βdt = [zero(t) for i in 1:max]
+    c_prime = Vector{typeof(t)}(undef, max) #time stepping
+    for i in 1 : max
+        c_prime[i] = zero(t)
+    end
 
     dw1 = zero(u)
     ubuff = zero(u)
-    dw2 = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
+    dw2 = [similar(u, Complex{eltype(u)}) for _ in 1 : (max - 1) ÷ 2]
     recursivefill!.(dw2, false)
-    cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
+    cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (max - 1) ÷ 2]
     recursivefill!.(cubuff, false)
-    dw = Vector{typeof(u)}(undef, num_stages - 1)
+    dw = [zero(u) for i in 1 : max]
 
-    cont = Vector{typeof(u)}(undef, num_stages)
-    for i in 1 : num_stages
-        cont[i] = zero(u)
-    end
+    cont = [zero(u) for i in 1:max]
 
-    derivatives = Matrix{typeof(u)}(undef, num_stages, num_stages)
-    for i in 1 : num_stages, j in 1 : num_stages
+    derivatives = Matrix{typeof(u)}(undef, max, max)
+    for i in 1 : max, j in 1 : max
         derivatives[i, j] = zero(u)
     end
 
     fsalfirst = zero(rate_prototype)
-    fw = Vector{typeof(rate_prototype)}(undef, num_stages)
-    ks = Vector{typeof(rate_prototype)}(undef, num_stages)
-    for i in 1: num_stages
-        ks[i] = fw[i] = zero(rate_prototype)
-    end
+    fw = [zero(rate_prototype) for i in 1 : max]
+    ks = [zero(rate_prototype) for i in 1 : max]
+
     k = ks[1]
 
     J, W1 = build_J_W(alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val(true))
@@ -622,7 +647,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         error("Non-concrete Jacobian not yet supported by AdaptiveRadau.")
     end
 
-    W2 = [similar(J, Complex{eltype(W1)}) for _ in 1 : (num_stages - 1) ÷ 2]
+    W2 = [similar(J, Complex{eltype(W1)}) for _ in 1 : (max - 1) ÷ 2]
     recursivefill!.(W2, false)
 
     du1 = zero(rate_prototype)
@@ -640,18 +665,18 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
 
     linsolve2 = [
         init(LinearProblem(W2[i], _vec(cubuff[i]); u0 = _vec(dw2[i])), alg.linsolve, alias_A = true, alias_b = true,
-            assumptions = LinearSolve.OperatorAssumptions(true)) for i in 1 : (num_stages - 1) ÷ 2]
+            assumptions = LinearSolve.OperatorAssumptions(true)) for i in 1 : (max - 1) ÷ 2]
 
     rtol = reltol isa Number ? reltol : zero(reltol)
     atol = reltol isa Number ? reltol : zero(reltol)
 
     AdaptiveRadauCache(u, uprev,
-        z, w, c_prime, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, 
+        z, w, c_prime, αdt, βdt, dw1, ubuff, dw2, cubuff, dw, cont, derivatives,
         du1, fsalfirst, ks, k, fw,
         J, W1, W2,
-        uf, tab, κ, one(uToltype), 10000, tmp,
+        uf, tabs, κ, one(uToltype), 10000, tmp,
         atmp, jac_config,
         linsolve1, linsolve2, rtol, atol, dt, dt,
-        Convergence, alg.step_limiter!)
+        Convergence, alg.step_limiter!, num_stages, 1, 0.0)
 end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
index 044257afad..2058c4fb15 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
@@ -90,7 +90,7 @@ function initialize!(integrator, cache::RadauIIA9Cache)
 end
 
 function initialize!(integrator, cache::AdaptiveRadauCache)
-    @unpack num_stages = cache.tab
+    @unpack num_stages = cache
     integrator.kshortsize = 2
     resize!(integrator.k, integrator.kshortsize)
     integrator.k[1] = integrator.fsalfirst
@@ -1032,7 +1032,7 @@ end
     @unpack dw1, ubuff, dw23, dw45, cubuff1, cubuff2 = cache
     @unpack k, k2, k3, k4, k5, fw1, fw2, fw3, fw4, fw5 = cache
     @unpack J, W1, W2, W3 = cache
-    @unpack tmp, tmp2, tmp3, tmp4, tmp5, tmp6, atmp, jac_config, linsolve1, linsolve2, linsolve3, rtol, atol, step_limiter! = cache
+    @unpack tmp, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, atmp, jac_config, linsolve1, linsolve2, linsolve3, rtol, atol, step_limiter! = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
     @unpack maxiters = alg
@@ -1087,30 +1087,30 @@ end
         c2′ = c2 * c5′
         c3′ = c3 * c5′
         c4′ = c4 * c5′
-        z1 = @.. c1′ * (cont1 +
+        @.. z1 = c1′ * (cont1 +
                                 (c1′-c4m1) * (cont2 +
                                                      (c1′ - c3m1) * (cont3 +
                                                                             (c1′ - c2m1) * (cont4 + (c1′ - c1m1) * cont5))))
-        z2 = @.. c2′ * (cont1 +
+        @.. z2 = c2′ * (cont1 +
                                 (c2′-c4m1) * (cont2 +
                                                      (c2′ - c3m1) * (cont3 +
                                                                             (c2′ - c2m1) * (cont4 + (c2′ - c1m1) * cont5))))
-        z3 = @.. c3′ * (cont1 +
+        @.. z3 = c3′ * (cont1 +
                                 (c3′-c4m1) * (cont2 +
                                                      (c3′ - c3m1) * (cont3 +
                                                                             (c3′ - c2m1) * (cont4 + (c3′ - c1m1) * cont5))))
-        z4 = @.. c4′ * (cont1 +
+        @.. z4 = c4′ * (cont1 +
                                 (c4′-c4m1) * (cont2 +
                                                      (c4′ - c3m1) * (cont3 +
                                                                             (c4′ - c2m1) * (cont4 + (c4′ - c1m1) * cont5))))
-        z5 = @.. c5′ * (cont1 +
+        @.. z5 = c5′ * (cont1 +
                                 (c5′-c4m1) * (cont2 +
                                                      (c5′ - c3m1) * (cont3 + (c5′ - c2m1) * (cont4 + (c5′ - c1m1) * cont5))))
-        w1 = @.. broadcast=false TI11*z1+TI12*z2+TI13*z3+TI14*z4+TI15*z5
-        w2 = @.. broadcast=false TI21*z1+TI22*z2+TI23*z3+TI24*z4+TI25*z5
-        w3 = @.. broadcast=false TI31*z1+TI32*z2+TI33*z3+TI34*z4+TI35*z5
-        w4 = @.. broadcast=false TI41*z1+TI42*z2+TI43*z3+TI44*z4+TI45*z5
-        w5 = @.. broadcast=false TI51*z1+TI52*z2+TI53*z3+TI54*z4+TI55*z5
+        @.. w1 = TI11*z1+TI12*z2+TI13*z3+TI14*z4+TI15*z5
+        @.. w2 = TI21*z1+TI22*z2+TI23*z3+TI24*z4+TI25*z5
+        @.. w3 = TI31*z1+TI32*z2+TI33*z3+TI34*z4+TI35*z5
+        @.. w4 = TI41*z1+TI42*z2+TI43*z3+TI44*z4+TI45*z5
+        @.. w5 = TI51*z1+TI52*z2+TI53*z3+TI54*z4+TI55*z5
     end
 
     # Newton iteration
@@ -1328,21 +1328,21 @@ end
     if integrator.EEst <= oneunit(integrator.EEst)
         cache.dtprev = dt
         if alg.extrapolant != :constant
-            cache.cont1 = @.. (z4 - z5) / c4m1 # first derivative on [c4, 1]
-            tmp1 = @.. (z3 - z4) / c3mc4 # first derivative on [c3, c4]
-            cache.cont2 = @.. (tmp1 - cache.cont1) / c3m1 # second derivative on [c3, 1]
-            tmp2 = @.. (z2 - z3) / c2mc3 # first derivative on [c2, c3]
-            tmp3 = @.. (tmp2 - tmp1) / c2mc4 # second derivative on [c2, c4]
-            cache.cont3 = @.. (tmp3 - cache.cont2) / c2m1 # third derivative on [c2, 1]
-            tmp4 = @.. (z1 - z2) / c1mc2 # first derivative on [c1, c2]
-            tmp5 = @.. (tmp4 - tmp2) / c1mc3 # second derivative on [c1, c3]
-            tmp6 = @.. (tmp5 - tmp3) / c1mc4 # third derivative on [c1, c4]
-            cache.cont4 = @.. (tmp6 - cache.cont3) / c1m1 #fourth derivative on [c1, 1]
-            tmp7 = @.. z1 / c1 #first derivative on [0, c1]
-            tmp8 = @.. (tmp4 - tmp7) / c2 #second derivative on [0, c2]
-            tmp9 = @.. (tmp5 - tmp8) / c3 #third derivative on [0, c3]
-            tmp10 = @.. (tmp6 - tmp9) / c4 #fourth derivative on [0,c4]
-            cache.cont5 = @.. cache.cont4 - tmp10 #fifth derivative on [0,1]
+            @.. cache.cont1 = (z4 - z5) / c4m1 # first derivative on [c4, 1]
+            @.. tmp = (z3 - z4) / c3mc4 # first derivative on [c3, c4]
+            @.. cache.cont2 = (tmp - cache.cont1) / c3m1 # second derivative on [c3, 1]
+            @.. tmp2 = (z2 - z3) / c2mc3 # first derivative on [c2, c3]
+            @.. tmp3 = (tmp2 - tmp) / c2mc4 # second derivative on [c2, c4]
+            @.. cache.cont3 = (tmp3 - cache.cont2) / c2m1 # third derivative on [c2, 1]
+            @.. tmp4 = (z1 - z2) / c1mc2 # first derivative on [c1, c2]
+            @.. tmp5 = (tmp4 - tmp2) / c1mc3 # second derivative on [c1, c3]
+            @.. tmp6 = (tmp5 - tmp3) / c1mc4 # third derivative on [c1, c4]
+            @.. cache.cont4 = (tmp6 - cache.cont3) / c1m1 #fourth derivative on [c1, 1]
+            @.. tmp7 = z1 / c1 #first derivative on [0, c1]
+            @.. tmp8 = (tmp4 - tmp7) / c2 #second derivative on [0, c2]
+            @.. tmp9 = (tmp5 - tmp8) / c3 #third derivative on [0, c3]
+            @.. tmp10 = (tmp6 - tmp9) / c4 #fourth derivative on [0,c4]
+            @.. cache.cont5 = cache.cont4 - tmp10 #fifth derivative on [0,1]
         end
     end
 
@@ -1354,7 +1354,9 @@ end
 @muladd function perform_step!(integrator, cache::AdaptiveRadauConstantCache,
     repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab
+    @unpack tabs, num_stages = cache
+    tab = tabs[(num_stages - 1) ÷ 2]
+    @unpack T, TI, γ, α, β, c, e = tab
     @unpack κ, cont = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
@@ -1398,16 +1400,16 @@ end
             cache.cont[i] = @.. map(zero, u)
         end
     else
-        c_prime = Vector{typeof(u)}(undef, num_stages) #time stepping
-        c_prime[num_stages] = @.. dt / cache.dtprev
+        c_prime = Vector{typeof(dt)}(undef, num_stages) #time stepping
+        c_prime[num_stages] = dt / cache.dtprev
         for i in 1 : num_stages - 1
-            c_prime[i] = @.. c[i] * c_prime[num_stages]
+            c_prime[i] = c[i] * c_prime[num_stages]
         end
         for i in 1 : num_stages # collocation polynomial
-            z[i] = @.. cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1]
+            z[i] = @.. cont[num_stages - 1] + cont[num_stages] * (c_prime[i] - c[1] + 1)
             j = num_stages - 2
             while j > 0
-                z[i] = @.. z[i] * (c_prime[i] - c[num_stages - j] + 1) + cont[j]
+                z[i] = @.. cont[j] + z[i] * (c_prime[i] - c[num_stages - j] + 1)
                 j = j - 1
             end
             z[i] = @.. z[i] * c_prime[i]
@@ -1431,45 +1433,45 @@ end
         integrator.stats.nnonliniter += 1
 
         # evaluate function
-        ff = Vector{typeof(u)}(undef, num_stages)
+        #ff = Vector{typeof(u)}(undef, num_stages)
         for i in 1 : num_stages
-            ff[i] = f(uprev + z[i], p, t + c[i] * dt)
+            z[i] = f(uprev + z[i], p, t + c[i] * dt)
         end
-        integrator.stats.nf += num_stages
+        OrdinaryDiffEqCore.increment_nf!(integrator.stats, num_stages)
 
         #fw = TI * ff
         fw = Vector{typeof(u)}(undef, num_stages)
         for i in 1 : num_stages
             fw[i] = @.. zero(u)
             for j in 1:num_stages
-                fw[i] = @.. fw[i] + TI[i,j] * ff[j]
+                fw[i] = @.. fw[i] + TI[i,j] * z[j]
             end
         end
 
-        Mw = Vector{typeof(u)}(undef, num_stages)
+        #Mw = Vector{typeof(u)}(undef, num_stages)
         if mass_matrix isa UniformScaling # `UniformScaling` doesn't play nicely with broadcast
             for i in 1 : num_stages
-                Mw[i] = @.. mass_matrix.λ * w[i] #scaling by eigenvalue
+                z[i] = @.. mass_matrix.λ * w[i] #scaling by eigenvalue
             end
         else
-            Mw = mass_matrix * w #standard multiplication
+            z = mass_matrix * w #standard multiplication
         end
 
         rhs = Vector{typeof(u)}(undef, num_stages)
-        rhs[1] = @.. fw[1] - γdt * Mw[1]
+        rhs[1] = @.. fw[1] - γdt * z[1]
         i = 2
         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]
+            rhs[i] = @.. fw[i] - αdt[i ÷ 2] * z[i] + βdt[i ÷ 2] * z[i + 1]
+            rhs[i + 1] = @.. fw[i + 1] - βdt[i ÷ 2] * z[i] - αdt[i ÷ 2] * z[i + 1]
             i += 2
         end
 
-        dw = Vector{typeof(u)}(undef, num_stages)
-        dw[1] = _reshape(LU1 \ _vec(rhs[1]), axes(u))
+        #dw = Vector{typeof(u)}(undef, num_stages)
+        z[1] = _reshape(LU1 \ _vec(rhs[1]), axes(u))
         for i in 2 :(num_stages + 1) ÷ 2
             tmp = _reshape(LU2[i - 1] \ _vec(@.. rhs[2 * i - 2] + rhs[2 * i - 1] * im), axes(u))
-            dw[2 * i - 2] = @.. real(tmp)
-            dw[2 * i - 1] = @.. imag(tmp)
+            z[2 * i - 2] = @.. real(tmp)
+            z[2 * i - 1] = @.. imag(tmp)
         end
         integrator.stats.nsolve +=(num_stages + 1) ÷ 2
 
@@ -1477,7 +1479,7 @@ end
         iter > 1 && (ndwprev = ndw)
         ndw = 0.0
         for i in 1 : num_stages
-            ndw += internalnorm(calculate_residuals(dw[i], uprev, u, atol, rtol, internalnorm, t), t)
+            ndw += internalnorm(calculate_residuals(z[i], uprev, u, atol, rtol, internalnorm, t), t)
         end
 
         # check divergence (not in initial step)
@@ -1493,7 +1495,7 @@ end
         end
         
         for i in 1 : num_stages
-            w[i] = @.. w[i] - dw[i]
+            w[i] = @.. w[i] - z[i]
         end
 
         # transform `w` to `z`
@@ -1534,10 +1536,13 @@ end
     u = @.. uprev + z[num_stages]
     
     if adaptive
-        edt = e ./ dt
-        tmp = dot(edt, z)
+        tmp = 0
+        for i in 1 : num_stages
+            tmp = @.. tmp + e[i]/dt * z[i]
+        end
         mass_matrix != I && (tmp = mass_matrix * tmp)
-        utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst+tmp
+        #utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst + tmp
+        utilde = @.. broadcast=false integrator.fsalfirst + tmp
         if alg.smooth_est
             utilde = _reshape(LU1 \ _vec(utilde), axes(u))
             integrator.stats.nsolve += 1
@@ -1548,8 +1553,9 @@ end
         if !(integrator.EEst < oneunit(integrator.EEst)) && integrator.iter == 1 ||
         integrator.u_modified
             f0 = f(uprev .+ utilde, p, t)
-            integrator.stats.nf += 1
-            utilde = @.. broadcast=false 1 / γ * dt * f0 + tmp
+            OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+            #utilde = @.. broadcast=false 1 / γ * dt * f0 + tmp
+            utilde = @.. broadcast=false f0+tmp
             if alg.smooth_est
                 utilde = _reshape(LU1 \ _vec(utilde), axes(u))
                 integrator.stats.nsolve += 1
@@ -1589,8 +1595,10 @@ 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, num_stages = cache.tab
-    @unpack κ, cont, derivatives, z, w, c_prime = cache
+    @unpack num_stages, tabs = cache
+    tab = tabs[(num_stages - 1) ÷ 2]
+    @unpack T, TI, γ, α, β, c, e = tab
+    @unpack κ, cont, derivatives, z, w, c_prime, αdt, βdt= cache
     @unpack dw1, ubuff, dw2, cubuff, dw = cache
     @unpack ks, k, fw, J, W1, W2 = cache
     @unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter! = cache
@@ -1600,13 +1608,18 @@ end
     mass_matrix = integrator.f.mass_matrix
 
     # precalculations
-    γdt, αdt, βdt = γ / dt, α ./ dt, β ./ dt
+    γdt = γ / dt
+    for i in 1 : (num_stages - 1) ÷ 2
+        αdt[i] = α[i]/dt
+        βdt[i] = β[i]/dt
+    end
+
     (new_jac = do_newJ(integrator, alg, cache, repeat_step)) &&
         (calc_J!(J, integrator, cache); cache.W_γdt = dt)
     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 :(num_stages - 1) ÷ 2
+            for i in 1 : (num_stages - 1) ÷ 2
                 W2[i][II] = -(αdt[i] + βdt[i] * im) * mass_matrix[Tuple(II)...] + J[II]
             end
         end
@@ -1630,7 +1643,8 @@ end
             @.. z[i] = cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1]
             j = num_stages - 2
             while j > 0
-                @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) + cont[j]
+                @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) 
+                @.. z[i] += cont[j]
                 j = j - 1
             end
             @.. z[i] *= c_prime[i]
@@ -1659,13 +1673,13 @@ end
             @.. tmp = uprev + z[i]
             f(ks[i], tmp, p, t + c[i] * dt)
         end
-        integrator.stats.nf += num_stages
+        OrdinaryDiffEqCore.increment_nf!(integrator.stats, num_stages)
 
         #mul!(fw, TI, ks)
         for i in 1:num_stages
-            fw[i] = @.. zero(u)
+            @.. fw[i] = zero(u)
             for j in 1:num_stages
-                fw[i] = @.. fw[i] + TI[i,j] * ks[j]
+                @.. fw[i] += TI[i,j] * ks[j]
             end
         end
 
@@ -1686,35 +1700,27 @@ end
         @.. ubuff = fw[1] - γdt * Mw[1]
         needfactor = iter == 1 && new_W
 
-        linsolve1 = cache.linsolve1
         if needfactor
-            linres = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff), linu = _vec(dw1))
+            cache.linsolve1 = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff), linu = _vec(dw1)).cache
         else
-            linres = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), linu = _vec(dw1))
+            cache.linsolve1 = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), linu = _vec(dw1)).cache
         end
 
-        cache.linsolve1 = linres.cache
-
-        linres2 = Vector{Any}(undef,(num_stages - 1) ÷ 2)
-
         for i in 1 :(num_stages - 1) ÷ 2
             @.. cubuff[i]=complex(
             fw[2 * i] - αdt[i] * Mw[2 * i] + βdt[i] * Mw[2 * i + 1], fw[2 * i + 1] - βdt[i] * Mw[2 * i] - αdt[i] * Mw[2 * i + 1])
             if needfactor
-                linres2[i] = dolinsolve(integrator, linsolve2[i]; A = W2[i], b = _vec(cubuff[i]), linu = _vec(dw2[i]))
+                cache.linsolve2[i] = dolinsolve(integrator, linsolve2[i]; A = W2[i], b = _vec(cubuff[i]), linu = _vec(dw2[i])).cache
             else
-                linres2[i] = dolinsolve(integrator, linsolve2[i]; A = nothing, b = _vec(cubuff[i]), linu = _vec(dw2[i]))
+                cache.linsolve2[i] = dolinsolve(integrator, linsolve2[i]; A = nothing, b = _vec(cubuff[i]), linu = _vec(dw2[i])).cache
             end
-            cache.linsolve2[i] = linres2[i].cache
         end
 
         integrator.stats.nsolve += (num_stages + 1) / 2
 
         for i in 1 : (num_stages - 1) ÷ 2
-            dw[2 * i - 1] = z[2 * i - 1]
-            dw[2 * i] = z[2 * i]
-            dw[2 * i - 1] = real(dw2[i])
-            dw[2 * i] = imag(dw2[i])
+            @.. dw[2 * i - 1] = real(dw2[i])
+            @.. dw[2 * i] = imag(dw2[i])
         end
 
         # compute norm of residuals
@@ -1746,15 +1752,15 @@ end
         # transform `w` to `z`
         #mul!(z, T, w)
         for i in 1:num_stages - 1
-            z[i] = @.. zero(u)
+            @.. z[i] = zero(u)
             for j in 1:num_stages
-                z[i] = @.. z[i] + T[i,j] * w[j]
+                @.. z[i] += T[i,j] * w[j]
             end
         end
-        z[num_stages] = @.. T[num_stages, 1] * w[1]
+        @.. z[num_stages] = T[num_stages, 1] * w[1]
         i = 2
         while i < num_stages
-            z[num_stages] = @.. z[num_stages] + w[i]
+            @.. z[num_stages] += w[i]
             i += 2
         end
 
@@ -1782,16 +1788,18 @@ end
     step_limiter!(u, integrator, p, t + dt)
     
     if adaptive
-        utilde = w2
-        edt = e./dt
-        @.. tmp = dot(edt, z) + 1 / γ * dt * fsalfirst
-        mass_matrix != I && (mul!(w1, mass_matrix, tmp); copyto!(tmp, w1))
-        @.. ubuff=integrator.fsalfirst + tmp
+        utilde = w[2]
+        @.. tmp = 0
+        for i in 1 : num_stages
+            @.. tmp += e[i]/dt * z[i]
+        end
+        mass_matrix != I && (mul!(w[1], mass_matrix, tmp); copyto!(tmp, w[1]))
+        #@.. ubuff=1 / γ * dt * integrator.fsalfirst + tmp
+        @.. broadcast=false ubuff=integrator.fsalfirst + tmp
 
         if alg.smooth_est
-            linres1 = dolinsolve(integrator, linres1.cache; b = _vec(ubuff),
-                linu = _vec(utilde))
-            cache.linsolve1 = linres1.cache
+            cache.linsolve1 = dolinsolve(integrator, linsolve1; b = _vec(ubuff),
+            linu = _vec(utilde)).cache
             integrator.stats.nsolve += 1
         end
 
@@ -1804,13 +1812,13 @@ end
            integrator.u_modified
             @.. broadcast=false utilde=uprev + utilde
             f(fsallast, utilde, p, t)
-            integrator.stats.nf += 1
+            OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+            #@.. broadcast=false ubuff = 1 / γ * dt * fsallast + tmp
             @.. broadcast=false ubuff=fsallast + tmp
 
             if alg.smooth_est
-                linres1 = dolinsolve(integrator, linres1.cache; b = _vec(ubuff),
-                    linu = _vec(utilde))
-                cache.linsolve1 = linres1.cache
+                cache.linsolve1 = dolinsolve(integrator, linsolve1; b = _vec(ubuff),
+                linu = _vec(utilde)).cache
                 integrator.stats.nsolve += 1
             end
 
@@ -1833,7 +1841,7 @@ end
                 end
             end
             for i in 1 : num_stages
-                cache.cont[i] = derivatives[i, num_stages]
+                @.. cache.cont[i] = derivatives[i, num_stages]
             end
         end
     end
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index 750a9a0b34..7191eabdb2 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -43,7 +43,7 @@ struct RadauIIA5Tableau{T, T2}
     T22::T
     T23::T
     T31::T
-    #T32::T 
+    #T32::T
     #T33::T
     TI11::T
     TI12::T
@@ -111,141 +111,6 @@ function RadauIIA5Tableau(T, T2)
         e1, e2, e3)
 end
 
-struct RadauIIATableau{T1, T2}
-    T::Matrix{T1}
-    TI::Matrix{T1}
-    c::Vector{T2}
-    γ::T1
-    α::Vector{T1}
-    β::Vector{T1}
-    e::Vector{T1}
-    num_stages::Int
-end
-
-function BigRadauIIA5Tableau(T1, T2)
-    γ = convert(T1, big"3.63783425274449573220841851357777579794593608687391153215117488565841871456727143375130115708511223004183651123208497057248238260532214672028700625775335843")
-    α = Vector{T1}(undef, 1)
-    β = Vector{T1}(undef, 1)
-    α[1] = big"2.68108287362775213389579074321111210102703195656304423392441255717079064271636428312434942145744388497908174438395751471375880869733892663985649687112332242"
-    β[1] = big"3.05043019924741056942637762478756790444070419917947659226291744751211727051786694870515117615266028855554735929171362769761399150862332538376382934625577549"
-
-    c = Vector{T2}(undef, 3)
-    c[1] = big"0.155051025721682190180271592529410860803405251934332987156730743274903962254268497346014056689535976518140539877338581087514113454016224265837421604876272084"
-    c[2] = big"0.644948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023481859460122661418912485886545983775734162578395123729143"
-    c[3] = big"1"
-
-    e = Vector{T1}(undef, 3)
-    e[1] = big"-0.428298294115368098113417591057340987284723986878723769598436582629558514867595819404541110575367601847354683220540647741034052880983125451920841851066713815"
-    e[2] = big"0.245039074384916438547779982042524963814308131639809054920681599475247366387530452311400517264433163448821319862243292031101333835037542080594323438762615605"
-    e[3] = big"-0.0916296098652259003910911359018870983677439465358993542342383692664256286871842771687905889462502859518430848371143253189423645209875225522045944109790661043"
-
-    TI = Matrix{T1}(undef, 3, 3)
-    TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
-    TI[1, 2] = big"0.339199251815809869542824974053410987511771566126056902312311333553438988409693737874718833892037643701271502187763370262948704203562215007824701228014200056"
-    TI[1, 3] = big"0.541770539935874871186523033492089631898841317849243944095021379289933921771713116368931784890546144473788347538203807242114936998948954098533375649163016612"
-    TI[2, 1] = big"-4.17871859155190472734646265851205623000038388214686525896709481539843195209360778128456932548583273459040707932166364293012713818843609182148794380267482041"
-    TI[2, 2] = big"-0.327682820761062387082533272429616234245791838308340887801415258608836530255609335712523838667242449344879454518796849992049787172023800373390124427898159896"
-    TI[2, 3] = big"0.476623554500550451960069084091012497939942928625055897109833707684876604712862299049343675491204859381277636585708398915065951363736337328178192801074535132"
-    TI[3, 1] = big"-0.502872634945786875951247343139544292859248429570937886791036339034110181540695221500843782634464164585836226038438397328726973424362168221527501738985822875"
-    TI[3, 2] = big"2.57192694985560542918678535360167505469448742842178326395573566888176471664393761903447163100353067504020263109067033226021288356347565113471227052083596358"
-    TI[3, 3] = big"-0.596039204828224924968821911099302403289857517521591823052174732952989090998130905722763344484798508456930766594977798579939415052669401095404149917833710127"
-
-    T = Matrix{T1}(undef, 3, 3)
-    T[1, 1] = big"0.091232394870892942791548135249436196118684699372210280712184363514099824021240149574725365814781580305065489937969163922775110463056339192206701819661425186"
-    T[1, 2] = big"-0.141255295020954208427990383807797309409263248498594798844289981408804297900674604638610419147468875667691398225003133444988034605081071965848437945842767211"
-    T[1 ,3] = big"-0.0300291941051474244918611170890538666683842974606300802563717702200388818691214144173874588956764952224874407424115249418136547481236684478531215095064078994"
-    T[2, 1] = big"0.241717932707107018957474779310148232884879540532595279746187345714229132659465207414913313803429072060469564350914390845001169448350326344874859416624577348"
-    T[2, 2] = big"0.204129352293799931995990810298338174086540402523315938937516234649384944528706774788799548853122282827246947911905379230680096946800308693162079538975632443"
-    T[2, 3] = big"0.382942112757261937795438233599873210357792575012007744255205163027042915338009760005422153613194350161760232119048691964499888989151661861236831969497483828"
-    T[3, 1] = big"0.966048182615092936190567080794590794996748754810883844283183333914131408744555961195911605614405476210484499875001737558078500322423463946527349731087504518"
-    T[3, 2] = big"1.0"
-    T[3, 3] = big"0.0"
-    RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, e, 3)
-end
-
-function BigRadauIIA9Tableau(T1, T2)
-    γ = convert(T1, big"6.28670475172927664517315334186940904959068186655567041187229167532923622489525703260842273089261139845280626287956099768662193453067483410165932355981736786")
-    α = Vector{T1}(undef, 2)
-    β = Vector{T1}(undef, 2)
-    α[1] = big"3.65569432546357225824320796009543385435699888857815445045567025741630720509235614026228963385258117304229337679733945535812317372403535763551850772878775217"
-    α[2] = big"5.70095329867178941917021536896986162084766017814401034360818390491907468246001534343349900070111312773130349176288004579856585901062722531365183049130382405"
-    β[1] = big"6.5437368993600772940210715093936863183637851728134458820202187133882261290012752452972782843700946890488789462524897903624959996932392239962196563965573345"
-    β[2] = big"3.21026560030854988842501065297211721232153653493981008029923647488964744732168461657389754087826565709085773529539707072244537983491480773006949966789260925"
-
-    c = Vector{T2}(undef, 5)
-    c[1] = big"0.0571041961145176821931211925541156212350779455987501643278082929309346782020731645861138168198427368635148018903413155731609901559772929443100370500757072557"
-    c[2] = big"0.276843013638123827680045997685625141110889169695030468349442048831121339683708036772541528564051130879197377136636984534220758899839905855114024309075271826"
-    c[3] = big"0.583590432368916820056697668662917248693432639896771640176293841831747501961831012005632277467456299345321045569611992496682381919275766424103024358378365496"
-    c[4] = big"0.8602401356562194478479129188751197667383780225872255049242335941839742579301655644134901549264276106897445531811874851737136468026848125542506920602484255"
-    c[5] = big"1.0"
-
-    e = Vector{T1}(undef, 5)
-    e[1] = big"-0.239909571163200476817707991076962793618683493830916562279975225042872448414819259070978815977101189851237591634144816820425592764432710089981892023441549743"
-    e[2] = big"0.125293484229223300606887443525929336197638450194973243323835481816625682669684634271160851996902445722310139152840852195000603355301430153561918341655137084"
-    e[3] = big"-0.0688288849083782089370741375422279772873399871312158026536514369967825732836040693366396751988019623495452730460577336089848105495733304937016477629621990433"
-    e[4] = big"0.0372433600301293198267284480468961585078575499935477902539140092558239369583143189611737274891328175087350382605569395937945872776839987254373869550943049195"
-    e[5] = big"-0.012863950751139890646895902730137465239579845437088427263654957605488133042638147254426913683751171160691603649073170415735165315443538347036196755548109703"
-
-    TI = Matrix{T1}(undef, 5, 5)
-    TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
-    TI[1, 2] = big"13.8651078562714131651762946846279728486098595017962436746405940971751244384714668104145151259298432908422191238542910724677205181071665482818120092330632702"
-    TI[1 ,3] = big"3.48000277479518556182840016971955819123081637245954095062693470191383865922357339844125383481645392882289968250993872221445874555610460465838129969397069557"
-    TI[1, 4] = big"-1.03200879782526342277108071214631493513824682491749273908106331923801396656058254294323988505859654767877050109789490714699847664805679842903430004696170252"
-    TI[1, 5] = big"0.804303045073989917475330383606196086089578671788707543063308602519859970319818304759856653218877415405946945572102875643297890954688508528143272905631829894"
-    TI[2, 1] = big"5.34418643783491159889531030409736033885455686563071401172022718575590068536629704134603404624953791012861634674294690788961703408019660066685859393456498931"
-    TI[2, 2] = big"4.59361556775916100445407449817656238428260055301676371438973411021009514435572975394999086474831271997070798032181411537895658457000537727156665947774751386"
-    TI[2, 3] = big"-3.03636032345942429864615756872018980250277648141683630832856906288036929718223473102394179699607901856890769270810252103326382063852039607285826867723587514"
-    TI[2, 4] = big"1.05066019023145886385983615715299311307615150447133905233370933194949591737765763708886464382722316727972166443876395823044171403663404254906698768838255919"
-    TI[2, 5] = big"-0.272778611864296270538614649997366804891835224042737605275699398413256470423268908248569612750117948720141667949532252500428432062582365619208502333677907158"
-    TI[3, 1] = big"3.74805980743980486005103450189256983678052751095791526209741655305580351377124372457009580386663275146166007984852101733055495783906881063060757645038080343"
-    TI[3, 2] = big"-3.98496573634388466725226385805351110838575115293851360514636734529255361185420464416807882769853298186283398369873418552760618971047757002216338511286260041"
-    TI[3, 3] = big"-1.04441564160801879294224732309562532189841624726401645191058551173485917137499204844819781779667611903670073971659834929382224472890100209497741235960707456"
-    TI[3, 4] = big"1.18409856813794848723102038838340482030291345603197522521517834943166421242518751666675199211369552058487095283489346390066317584532997854692445653563909898"
-    TI[3, 5] = big"-0.449917770156780368898811918314095435942113881883174152777026977062686286863549565130412864190301081537983106397709991028107600781961279985605930655683680139"
-    TI[4, 1] = big"-33.0418802135190000080614469426109507742858088371383868670878639187564531424382858814386742148456699143328462132296293097447566408853495288807407929988004676"
-    TI[4, 2] = big"-17.3769534790635670194549806058987105852733409102703844354448800193942184746909147697382687117638715195698950138089979798321855885541817752366521518811413713"
-    TI[4, 3] = big"-0.172129063254005561151528806427751383749451500597823574207174433146207178559871803504021077429693091164540897873472803934375603405253541639437370184767553293"
-    TI[4, 4] = big"-0.0991697779825426425881662214017368584726354746776989845479783944003623924121748016326495070834800297497011104846871751430208559227945252758721362340763610828"
-    TI[4, 5] = big"0.531228115838306667184911422606024795426589562580669892779793097035561488973256023529352389498509937781553683467106048413485632583844632286562240161995145055"
-    TI[5, 1] = big"-8.61144397987529197770008251257034851950485933115010902789613925540488896812417081206983938638600226846804467531843522104806738090683710882069500386691775154"
-    TI[5, 2] = big"9.69999140952880823133589405342003266497120753048627084327055311528684684237122654108691149692242002085965723391934376924400492239317026460192827344970015484"
-    TI[5, 3] = big"1.91472863969687428485137560339172471528025297511003983469957355306260543484472462223194401768126877615795915146192537091374017807611943419264038682143890747"
-    TI[5, 4] = big"2.41869200608494002642656343408298350771199306961305597858229870375990977712805399625496435641846363295393762353024017195444763964531237381728801981679934304"
-    TI[5, 5] = big"-1.0474634879353374186944329992117360176590042540536055452919974336199826846201614544718272622833822842591012529895091659029452542118642301415759073410771819"
-
-    T = Matrix{T1}(undef, 5, 5)
-    T[1, 1] = big"0.0125175862205010458901356760368001462557655123420858705973577952199246108029451084239310924615007306721702298573083400752464277227557045438770401832498107968"
-    T[1, 2] = big"-0.0102420478179088270700863300668590125015813934827825923708366359399562125950804289592272678367034071306578383319296130180550178248531589487456925441921649293"
-    T[1 ,3] = big"0.0476738772902957238631839478592069782970238490568258436986723993118380988311441474394156362952631834786373081794857384127209450988829840886524135970873769918"
-    T[1, 4] = big"-0.0114785152552295147079415554121555049385506204591245712490409384029671974157542450636658532835395855844059342442518520033304129991000509527123870917346017759"
-    T[1, 5] = big"-0.0140198588928754102810778942934959307831026572823203692568448424056201483917805257790275956734469193171917730378117501915144713896813544630288006687542182225"
-    T[2, 1] = big"0.00149167015189538242900444775236282223594625052328927847572623038484966999313257893341818287477809424303168766872838075463220122499449382436194198620498144296"
-    T[2, 2] = big"0.050172864517371058162991380262646513853120568882725793734131676894272706020317186004736779675826101816279321643304301437029912742375638648226701787880031719"
-    T[2, 3] = big"-0.0943318191816114369806569003363724471884924328367212069321438749304281980331334016578193750445513659941246363262225907407726099492713722343006925656625258579"
-    T[2, 4] = big"-0.00766883074918016288515687679203608074116106558796378201472238095295554979920808799930579174190884587422912077296093093698836937450535804218413704866981728518"
-    T[2, 5] = big"0.024708578426518526812525205377780382655366504554979744093019395818934704623702078004474076773426928900579988063099593288435684744957695210778788200213260272"
-    T[3, 1] = big"0.072981876388087148622657299703669587832652508881663282287850495621401398441897288250625556038835308015912409648841893161563884759791665776933761278383553608"
-    T[3, 2] = big"-0.230539534043417946721421862180000422679228296568599014834226319726930529322581417981617275287468418138394077987361681288909676234537699721082090802790143303"
-    T[3, 3] = big"0.102703045380125899792210456947141185148813233939327773583525878521508211077874610560448598369259541346968946573971195783374996178436435357335759255990489434"
-    T[3, 4] = big"0.0193984639988289509112232896408330872285824216708905773930244363652651247181543158008567311548336143384128605013911312875018664026371225431993252265128272262"
-    T[3, 5] = big"0.0818003537037511708363908122287572533071340646031113975848869261019231448226334426630664318901554550460201409321555775999869184033436795623062614812355590017"
-    T[4, 1] = big"0.380091440003568104126439184355215575526619121262253024859378518379910007234696730891540745160675744992320824590679292148769326540463161583672773762554445506"
-    T[4, 2] = big"0.377893902248861249543862293745933995234687511602719536459666284734445918178134851270924212812363352965391508894581698067329905034837778770261095647458874628"
-    T[4, 3] = big"0.466744130332494359289559582964906703283968612669234331018678042733321473730897217606173184300477207393539851157929838664168404778962779344509707214938022808"
-    T[4, 4] = big"0.40760117128019906662166237021895987274626181127101561893104166874567447589187790736078997321464949349935802836110699884016973990503134772720646054039223561"
-    T[4, 5] = big"0.199682427886802525936540566022390695167018315867216115995143539347975271751460199398235415129329119718414206048034051939441434136353381864781262773401023899"
-    T[5, 1] = big"0.921978973681210488488254647415676321266345412943047462855852351388222898143904205962703147998267738964059170225806964893009202287585991334322032058414768529"
-    T[5, 2] = big"1.0"
-    T[5, 3] = big"0.0"
-    T[5, 4] = big"1.0"
-    T[5, 5] = big"0.0"
-
-    RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, e, 5)
-end
-
-
 struct RadauIIA9Tableau{T, T2}
     T11::T
     T12::T
@@ -394,34 +259,156 @@ function RadauIIA9Tableau(T, T2)
         e1, e2, e3, e4, e5)
 end
 
-function BigRadauIIA13Tableau(T1, T2)
+struct RadauIIATableau{T1, T2}
+    T::Matrix{T1}
+    TI::Matrix{T1}
+    c::Vector{T2}
+    γ::T1
+    α::Vector{T1}
+    β::Vector{T1}
+    e::Vector{T1}
+end
+
+function RadauIIATableau5(T1, T2)
+    γ = convert(T1, big"3.63783425274449573220841851357777579794593608687391153215117488565841871456727143375130115708511223004183651123208497057248238260532214672028700625775335843")
+    α = T1[big"2.68108287362775213389579074321111210102703195656304423392441255717079064271636428312434942145744388497908174438395751471375880869733892663985649687112332242"]
+    β = T2[big"3.05043019924741056942637762478756790444070419917947659226291744751211727051786694870515117615266028855554735929171362769761399150862332538376382934625577549"]
+
+    c = T2[big"0.155051025721682190180271592529410860803405251934332987156730743274903962254268497346014056689535976518140539877338581087514113454016224265837421604876272084",
+           big"0.644948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023481859460122661418912485886545983775734162578395123729143",
+           1]
+
+    e = T1[big"-10.0488093998274155624603295076470799145872107881988969663429493235855742140670683952596720105774938812433874028620997746246706860729547671304601625528869782",
+           big"1.38214273316074889579366284098041324792054412153223029967628265691890754740040172859300534391082721457672073619543310795800401940628810046379349588622031217",
+           -1//3]
+
+    TI = Matrix{T1}(undef, 3, 3)
+    TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
+    TI[1, 2] = big"0.339199251815809869542824974053410987511771566126056902312311333553438988409693737874718833892037643701271502187763370262948704203562215007824701228014200056"
+    TI[1, 3] = big"0.541770539935874871186523033492089631898841317849243944095021379289933921771713116368931784890546144473788347538203807242114936998948954098533375649163016612"
+    TI[2, 1] = big"-4.17871859155190472734646265851205623000038388214686525896709481539843195209360778128456932548583273459040707932166364293012713818843609182148794380267482041"
+    TI[2, 2] = big"-0.327682820761062387082533272429616234245791838308340887801415258608836530255609335712523838667242449344879454518796849992049787172023800373390124427898159896"
+    TI[2, 3] = big"0.476623554500550451960069084091012497939942928625055897109833707684876604712862299049343675491204859381277636585708398915065951363736337328178192801074535132"
+    TI[3, 1] = big"-0.502872634945786875951247343139544292859248429570937886791036339034110181540695221500843782634464164585836226038438397328726973424362168221527501738985822875"
+    TI[3, 2] = big"2.57192694985560542918678535360167505469448742842178326395573566888176471664393761903447163100353067504020263109067033226021288356347565113471227052083596358"
+    TI[3, 3] = big"-0.596039204828224924968821911099302403289857517521591823052174732952989090998130905722763344484798508456930766594977798579939415052669401095404149917833710127"
+
+    T = Matrix{T1}(undef, 3, 3)
+    T[1, 1] = big"0.091232394870892942791548135249436196118684699372210280712184363514099824021240149574725365814781580305065489937969163922775110463056339192206701819661425186"
+    T[1, 2] = big"-0.141255295020954208427990383807797309409263248498594798844289981408804297900674604638610419147468875667691398225003133444988034605081071965848437945842767211"
+    T[1 ,3] = big"-0.0300291941051474244918611170890538666683842974606300802563717702200388818691214144173874588956764952224874407424115249418136547481236684478531215095064078994"
+    T[2, 1] = big"0.241717932707107018957474779310148232884879540532595279746187345714229132659465207414913313803429072060469564350914390845001169448350326344874859416624577348"
+    T[2, 2] = big"0.204129352293799931995990810298338174086540402523315938937516234649384944528706774788799548853122282827246947911905379230680096946800308693162079538975632443"
+    T[2, 3] = big"0.382942112757261937795438233599873210357792575012007744255205163027042915338009760005422153613194350161760232119048691964499888989151661861236831969497483828"
+    T[3, 1] = big"0.966048182615092936190567080794590794996748754810883844283183333914131408744555961195911605614405476210484499875001737558078500322423463946527349731087504518"
+    T[3, 2] = 1
+    T[3, 3] = 0
+    RadauIIATableau{T1, T2}(T, TI,
+    c, γ, α, β, e)
+end
+
+function RadauIIATableau9(T1, T2)
+    γ = convert(T1, big"6.28670475172927664517315334186940904959068186655567041187229167532923622489525703260842273089261139845280626287956099768662193453067483410165932355981736786")
+
+    α = T1[big"3.65569432546357225824320796009543385435699888857815445045567025741630720509235614026228963385258117304229337679733945535812317372403535763551850772878775217",
+           big"5.70095329867178941917021536896986162084766017814401034360818390491907468246001534343349900070111312773130349176288004579856585901062722531365183049130382405"]
+    β = T1[big"6.5437368993600772940210715093936863183637851728134458820202187133882261290012752452972782843700946890488789462524897903624959996932392239962196563965573345",
+           big"3.21026560030854988842501065297211721232153653493981008029923647488964744732168461657389754087826565709085773529539707072244537983491480773006949966789260925"]
+
+    c = T2[big"0.0571041961145176821931211925541156212350779455987501643278082929309346782020731645861138168198427368635148018903413155731609901559772929443100370500757072557",
+           big"0.276843013638123827680045997685625141110889169695030468349442048831121339683708036772541528564051130879197377136636984534220758899839905855114024309075271826",
+           big"0.583590432368916820056697668662917248693432639896771640176293841831747501961831012005632277467456299345321045569611992496682381919275766424103024358378365496",
+           big"0.8602401356562194478479129188751197667383780225872255049242335941839742579301655644134901549264276106897445531811874851737136468026848125542506920602484255",
+           1.0]
+
+    e = T1[big"-27.7809339440646373047872078172168798923674228687740760060378492475924178050505976287227228556471699142365371740120443650701118024100678675823465762727483305",
+           big"3.64147849804921315271165508774289722904088750334220956841022786858917594981395319605788667956024462601802006251583142928630101075351336314632135787805261686",
+           big"-1.25254772116911872049065249430114914889315244289570569309128740586057170336299694248256681515155624683225624015343224399700466177251702555220815764199263189",
+           big"0.592003167184542872566205223775131812219687808327572130718908784863813558599641375147402991238481535050773351649645179780815453429071529988233376036688329872",
+           -1//5]
+
+    TI = Matrix{T1}(undef, 5, 5)
+    TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
+    TI[1, 2] = big"13.8651078562714131651762946846279728486098595017962436746405940971751244384714668104145151259298432908422191238542910724677205181071665482818120092330632702"
+    TI[1 ,3] = big"3.48000277479518556182840016971955819123081637245954095062693470191383865922357339844125383481645392882289968250993872221445874555610460465838129969397069557"
+    TI[1, 4] = big"-1.03200879782526342277108071214631493513824682491749273908106331923801396656058254294323988505859654767877050109789490714699847664805679842903430004696170252"
+    TI[1, 5] = big"0.804303045073989917475330383606196086089578671788707543063308602519859970319818304759856653218877415405946945572102875643297890954688508528143272905631829894"
+    TI[2, 1] = big"5.34418643783491159889531030409736033885455686563071401172022718575590068536629704134603404624953791012861634674294690788961703408019660066685859393456498931"
+    TI[2, 2] = big"4.59361556775916100445407449817656238428260055301676371438973411021009514435572975394999086474831271997070798032181411537895658457000537727156665947774751386"
+    TI[2, 3] = big"-3.03636032345942429864615756872018980250277648141683630832856906288036929718223473102394179699607901856890769270810252103326382063852039607285826867723587514"
+    TI[2, 4] = big"1.05066019023145886385983615715299311307615150447133905233370933194949591737765763708886464382722316727972166443876395823044171403663404254906698768838255919"
+    TI[2, 5] = big"-0.272778611864296270538614649997366804891835224042737605275699398413256470423268908248569612750117948720141667949532252500428432062582365619208502333677907158"
+    TI[3, 1] = big"3.74805980743980486005103450189256983678052751095791526209741655305580351377124372457009580386663275146166007984852101733055495783906881063060757645038080343"
+    TI[3, 2] = big"-3.98496573634388466725226385805351110838575115293851360514636734529255361185420464416807882769853298186283398369873418552760618971047757002216338511286260041"
+    TI[3, 3] = big"-1.04441564160801879294224732309562532189841624726401645191058551173485917137499204844819781779667611903670073971659834929382224472890100209497741235960707456"
+    TI[3, 4] = big"1.18409856813794848723102038838340482030291345603197522521517834943166421242518751666675199211369552058487095283489346390066317584532997854692445653563909898"
+    TI[3, 5] = big"-0.449917770156780368898811918314095435942113881883174152777026977062686286863549565130412864190301081537983106397709991028107600781961279985605930655683680139"
+    TI[4, 1] = big"-33.0418802135190000080614469426109507742858088371383868670878639187564531424382858814386742148456699143328462132296293097447566408853495288807407929988004676"
+    TI[4, 2] = big"-17.3769534790635670194549806058987105852733409102703844354448800193942184746909147697382687117638715195698950138089979798321855885541817752366521518811413713"
+    TI[4, 3] = big"-0.172129063254005561151528806427751383749451500597823574207174433146207178559871803504021077429693091164540897873472803934375603405253541639437370184767553293"
+    TI[4, 4] = big"-0.0991697779825426425881662214017368584726354746776989845479783944003623924121748016326495070834800297497011104846871751430208559227945252758721362340763610828"
+    TI[4, 5] = big"0.531228115838306667184911422606024795426589562580669892779793097035561488973256023529352389498509937781553683467106048413485632583844632286562240161995145055"
+    TI[5, 1] = big"-8.61144397987529197770008251257034851950485933115010902789613925540488896812417081206983938638600226846804467531843522104806738090683710882069500386691775154"
+    TI[5, 2] = big"9.69999140952880823133589405342003266497120753048627084327055311528684684237122654108691149692242002085965723391934376924400492239317026460192827344970015484"
+    TI[5, 3] = big"1.91472863969687428485137560339172471528025297511003983469957355306260543484472462223194401768126877615795915146192537091374017807611943419264038682143890747"
+    TI[5, 4] = big"2.41869200608494002642656343408298350771199306961305597858229870375990977712805399625496435641846363295393762353024017195444763964531237381728801981679934304"
+    TI[5, 5] = big"-1.0474634879353374186944329992117360176590042540536055452919974336199826846201614544718272622833822842591012529895091659029452542118642301415759073410771819"
+
+    T = Matrix{T1}(undef, 5, 5)
+    T[1, 1] = big"0.0125175862205010458901356760368001462557655123420858705973577952199246108029451084239310924615007306721702298573083400752464277227557045438770401832498107968"
+    T[1, 2] = big"-0.0102420478179088270700863300668590125015813934827825923708366359399562125950804289592272678367034071306578383319296130180550178248531589487456925441921649293"
+    T[1 ,3] = big"0.0476738772902957238631839478592069782970238490568258436986723993118380988311441474394156362952631834786373081794857384127209450988829840886524135970873769918"
+    T[1, 4] = big"-0.0114785152552295147079415554121555049385506204591245712490409384029671974157542450636658532835395855844059342442518520033304129991000509527123870917346017759"
+    T[1, 5] = big"-0.0140198588928754102810778942934959307831026572823203692568448424056201483917805257790275956734469193171917730378117501915144713896813544630288006687542182225"
+    T[2, 1] = big"0.00149167015189538242900444775236282223594625052328927847572623038484966999313257893341818287477809424303168766872838075463220122499449382436194198620498144296"
+    T[2, 2] = big"0.050172864517371058162991380262646513853120568882725793734131676894272706020317186004736779675826101816279321643304301437029912742375638648226701787880031719"
+    T[2, 3] = big"-0.0943318191816114369806569003363724471884924328367212069321438749304281980331334016578193750445513659941246363262225907407726099492713722343006925656625258579"
+    T[2, 4] = big"-0.00766883074918016288515687679203608074116106558796378201472238095295554979920808799930579174190884587422912077296093093698836937450535804218413704866981728518"
+    T[2, 5] = big"0.024708578426518526812525205377780382655366504554979744093019395818934704623702078004474076773426928900579988063099593288435684744957695210778788200213260272"
+    T[3, 1] = big"0.072981876388087148622657299703669587832652508881663282287850495621401398441897288250625556038835308015912409648841893161563884759791665776933761278383553608"
+    T[3, 2] = big"-0.230539534043417946721421862180000422679228296568599014834226319726930529322581417981617275287468418138394077987361681288909676234537699721082090802790143303"
+    T[3, 3] = big"0.102703045380125899792210456947141185148813233939327773583525878521508211077874610560448598369259541346968946573971195783374996178436435357335759255990489434"
+    T[3, 4] = big"0.0193984639988289509112232896408330872285824216708905773930244363652651247181543158008567311548336143384128605013911312875018664026371225431993252265128272262"
+    T[3, 5] = big"0.0818003537037511708363908122287572533071340646031113975848869261019231448226334426630664318901554550460201409321555775999869184033436795623062614812355590017"
+    T[4, 1] = big"0.380091440003568104126439184355215575526619121262253024859378518379910007234696730891540745160675744992320824590679292148769326540463161583672773762554445506"
+    T[4, 2] = big"0.377893902248861249543862293745933995234687511602719536459666284734445918178134851270924212812363352965391508894581698067329905034837778770261095647458874628"
+    T[4, 3] = big"0.466744130332494359289559582964906703283968612669234331018678042733321473730897217606173184300477207393539851157929838664168404778962779344509707214938022808"
+    T[4, 4] = big"0.40760117128019906662166237021895987274626181127101561893104166874567447589187790736078997321464949349935802836110699884016973990503134772720646054039223561"
+    T[4, 5] = big"0.199682427886802525936540566022390695167018315867216115995143539347975271751460199398235415129329119718414206048034051939441434136353381864781262773401023899"
+    T[5, 1] = big"0.921978973681210488488254647415676321266345412943047462855852351388222898143904205962703147998267738964059170225806964893009202287585991334322032058414768529"
+    T[5, 2] = 1
+    T[5, 3] = 0
+    T[5, 4] = 1
+    T[5, 5] = 0
+
+    RadauIIATableau{T1, T2}(T, TI,
+    c, γ, α, β, e)
+end
+
+function RadauIIATableau13(T1, T2)
     γ = convert(T1, big"8.93683278840521633730209691330107970355008194433956657198414191417624969654351559268800871286734194720118970058657997472527299153742511021973612156231867783")
-    α = Vector{T1}(undef, 3)
-    β = Vector{T1}(undef, 3)
-    α[1] = big"4.37869356150680600252334919268856129165763746518197948235657247177701087073069907016715245914093899486193202405685779803686971216417800783050995450529391908"
-    α[2] = big"7.14105521918764010577498142571556804318193862372238812855726792587872300446315860222917039505087745633962330233504078264632719519730762016919715839787116038"
-    α[3] = big"8.51183482510294572305062092494533081338538293892584910309408864525614127653438453125967278937451257519784982331481143195416659686980181689042482631568989031"
-    β[1] = big"10.1696932837950116273183544188477298930096536824510223588525334625762336174947183926243705927725260475934351162622185429326813205432867247703480391692806137"
-    β[2] = big"6.62304592263927597062055811591186110468148199066707542227575094761515104946479159063603447729283770429494038962408904312215452856333028405675512985803584472"
-    β[3] = big"3.2810136243250588300359425270393915846791621918405321383787427650552081712406957205287551182809705166989352673500472974040971593568323836675590314648604458"
-
-    c = Vector{T2}(undef, 7)
-    c[1] = big"0.0293164271597848919720502769131649103737303925637149277869106839449360382416657787486309483651843695097273923248526200112627747993405898353736305552306269904"
-    c[2] = big"0.148078599668484291849976852495979212230248774808594461412594641801598386090878321806369397661747576057906341132861865305306667654594593138746653233717241913"
-    c[3] = big"0.336984690281154299097052972080775705197568750028473347122562968073691350512784060852409141173654482529393236826516171319486086447256539582972346127980810124"
-    c[4] = big"0.558671518771550132081393341805521940074368288965407825555747226117350122897421078323820052012282581935200398463518265914564420109615277886000739200777932339"
-    c[5] = big"0.769233862030054500916883360115645451837142143322295416166948169636548130573953285685200211542774367652885154701431860087378103033801830280742146083476036669"
-    c[6] = big"0.926945671319741114851873965819682011056172419542283252724467079656645202452528243814339480013587391545656707320049986592771178724621938506933715568048004783"
-    c[7] = big"1.0"
-
-    e = Vector{T1}(undef, 7)
-    e[1] = big"-0.171003707892600662399316289094713451418682228802554671094543075316220821099884263713032871578607576486535539488632407200766379971279557791263375813632287147"
-    e[2] = big"0.0934967172358652400317534533028674569308657324394316331629203486361371292312231403973668315582870547753526899857449840409175610009437530537219068836721686211"
-    e[3] = big"-0.0538908303114758775848180855518003793385120454908028879947132475960997563222416509683480350817114056356343433378213334826358824351525243402758389616051172681"
-    e[4] = big"0.03036786965048439581219923157368590250090409822952169396779792168990510618756404452728392288892998298088147691907128776816545685599760715439221674418662785"
-    e[5] = big"-0.0169792974425458224044481617230998766694942329759644144506911809530904808476835739189122151558601434810772520378036293579816345384682687602578758514350075723"
-    e[6] = big"0.00942688256820236884916415666439281573527695349338346787075550606528435808025071461023926432308932314282041885090975780812273515256740094297311541275151861292"
-    e[7] = big"-0.00331409873565629283448601269346047459594635696619041493081994712789731442072563377354737487903843138987115421498455722938021358621090485566426506726181005806"
+    α = T1[big"4.37869356150680600252334919268856129165763746518197948235657247177701087073069907016715245914093899486193202405685779803686971216417800783050995450529391908",
+           big"7.14105521918764010577498142571556804318193862372238812855726792587872300446315860222917039505087745633962330233504078264632719519730762016919715839787116038",
+           big"8.51183482510294572305062092494533081338538293892584910309408864525614127653438453125967278937451257519784982331481143195416659686980181689042482631568989031"]
+    β = T1[big"10.1696932837950116273183544188477298930096536824510223588525334625762336174947183926243705927725260475934351162622185429326813205432867247703480391692806137",
+           big"6.62304592263927597062055811591186110468148199066707542227575094761515104946479159063603447729283770429494038962408904312215452856333028405675512985803584472",
+           big"3.2810136243250588300359425270393915846791621918405321383787427650552081712406957205287551182809705166989352673500472974040971593568323836675590314648604458"]
+
+    c = T2[big"0.0293164271597848919720502769131649103737303925637149277869106839449360382416657787486309483651843695097273923248526200112627747993405898353736305552306269904",
+           big"0.148078599668484291849976852495979212230248774808594461412594641801598386090878321806369397661747576057906341132861865305306667654594593138746653233717241913",
+           big"0.336984690281154299097052972080775705197568750028473347122562968073691350512784060852409141173654482529393236826516171319486086447256539582972346127980810124",
+           big"0.558671518771550132081393341805521940074368288965407825555747226117350122897421078323820052012282581935200398463518265914564420109615277886000739200777932339",
+           big"0.769233862030054500916883360115645451837142143322295416166948169636548130573953285685200211542774367652885154701431860087378103033801830280742146083476036669",
+           big"0.926945671319741114851873965819682011056172419542283252724467079656645202452528243814339480013587391545656707320049986592771178724621938506933715568048004783",
+           1]
+
+    e = T1[big"-54.374436894128614514583710369683221528326818668136315170227649609831132483812209590903458627819914413600703287942266678601263304348350182019714004102122958",
+           big"7.00002400425918651204068363735192307633403862621907697222411411256593188888314837387690262103761082115674234000933589934965063951414231971808906314491204573",
+           big"-2.35566109198755719225604586775720723211163199654640573606711168106849118084357027539414093812951288166804790294091903523762277368547775099880612390898224076",
+           big"1.13228906610613438638449290827978318662460499026070073842612187085281352278780837966549916347601259689966925986653914736463076068138934273474363230390185871",
+           big"-0.646891326767358711867345222439989069591870662562921671446738173180691199552327090727940249497816198076028398716990245669520129053944261569921119452534594627",
+           big"0.387533385375352377424782057105854424214534853623007724234120623518712309680007346340280888076477218145510846867158055651267664035097674992751409157682864641",
+           -1//7]
 
     TI = Matrix{T1}(undef, 7, 7)
     TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676"
@@ -518,116 +505,66 @@ function BigRadauIIA13Tableau(T1, T2)
     T[6, 6] = big"0.521751945274765285294609453181807034209434470364856664246194441011327338299794536726049398636575212016960129143954076748520870645966241492966592488607495009"
     T[6, 7] = big"0.128071944635543894414114939510913357662538610722706228789484435811417614332529416514635125851744500940930818246509599119254761178392202724896572159336577251"
     T[7, 1] = big"0.881391578353818376313498879127399181693003124999819194603124949551827789004545406999549226388170693806014968936224161749923163222614460424501073405017519348"
-    T[7, 2] = big"1.0"
-    T[7, 3] = big"0.0"
-    T[7, 4] = big"1.0"
-    T[7, 5] = big"0.0"
-    T[7, 6] = big"1.0"
-    T[7, 7] = big"0.0"
-
-    RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, e, 7)
+    T[7, 2] = 1
+    T[7, 3] = 0
+    T[7, 4] = 1
+    T[7, 5] = 0
+    T[7, 6] = 1
+    T[7, 7] = 0
+
+    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e)
 end
 
-using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
-using Symbolics: variables, variable, unwrap
+import LinearAlgebra: eigen
+import FastGaussQuadrature: gaussradau
 
-function adaptiveRadauTableau(T1, 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, 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
-    radau_p = Polynomial{BigFloat}([tmp; tmp2])
-    for i in 1:(num_stages - 1)
-        radau_p = derivative(radau_p)
+function RadauIIATableau(T1, T2, num_stages::Int)
+    c = reverse!(1 .- gaussradau(num_stages, T1)[1])./2
+    if T1 == T2
+        c2 = c
+    else
+        c2 = reverse!(1 .- gaussradau(num_stages, T2)[1])./2
     end
-    c = real(roots(radau_p))
-    c[num_stages] = 1
-    c_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
+
+    c_powers = Matrix{T1}(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)
+        c_powers[i, 1] = 1
+        for j in 2 : num_stages
+            c_powers[i,j] = c[i]*c_powers[i,j-1]
         end
     end
-    inverse_c_powers = inv(c_powers)
-    c_q = Matrix{BigFloat}(undef, num_stages, num_stages)
+    c_q = Matrix{T1}(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
+            c_q[i,j] = c_powers[i,j] * c[i] / j
         end
     end
-    a = c_q * inverse_c_powers
-    a_inverse = inv(a)
-    b = Vector{BigFloat}(undef, num_stages)
-    for i in 1 : num_stages
-        b[i] = a[num_stages, i]
-    end
-    vals = eigvals(a_inverse)
-    γ = real(vals[num_stages])
-    α = Vector{BigFloat}(undef, floor(Int, num_stages/2))
-    β = Vector{BigFloat}(undef, floor(Int, num_stages/2))
-    index = 1
-    i = 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, num_stages)
-    i = 1
-    index = 2
-    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
+    a = c_q / c_powers
+
+    local eigval, eigvec;
+    try
+        eigval, eigvec = eigen(a)
+    catch
+        throw(ArgumentError("Solving ODEs with AdaptiveRadau with $T1 eltype and max_order >=17 requires loading GenericSchur.jl"))
     end
-    vecs[1] = real(eigvec[:, num_stages])
-    tmp3 = vcat(vecs)
-    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
+    # α, β, and γ come from eigvals(inv(a)) which are equal to inv.(eivals(a))
+    eigval .= inv.(eigval)
+    α = [real(eigval[i]) for i in 1:2:num_stages-1]
+    β = [imag(eigval[i]) for i in 1:2:num_stages-1]
+    γ = real(eigval[num_stages])
+
+    T = Matrix{T1}(undef, num_stages, num_stages)
+    @views for i in 2:2:num_stages
+        T[:, i] .= real.(eigvec[:, i] ./ eigvec[num_stages, i])
+        T[:, i + 1] .= imag.(eigvec[:, i] ./ eigvec[num_stages, i])
     end
+    @views T[:, 1] .= real.(eigvec[:, num_stages])
     TI = inv(T)
-
-    if (num_stages == 9)
-        e = Vector{BigFloat}(undef, 9)
-        e[1] = big"-0.133101731359431287515066981129913748644705107621439651956220186897253838380345034218235538734967567153163030284540660584311040323114847240173627907922903296"
-        e[2] = big"0.0754476228408557299650196603226967248368445025181771896522057250989188754588885465998346476425502117889420021664297319179240040109156780754680742172762707621"
-        e[3] = big"-0.0458369394236156144604575482137179697005739995740615341890112217655441769701945378217626766299683076189687755618065050383493055018324395934911567207485032988"
-        e[4] = big"0.0271430329153098694457979735602502142083095152399102869109830450899844979409229538982100527256348792152825816553434603418662939944133319974874915933773657075"
-        e[5] = big"-0.0156126300301219212217568535995825232086423550686814635293876744035364259647929167763641353639085929285192729729570945658304937255929114458885296622493040224"
-        e[6] = big"0.00890598154557403928205152521539967562877335780940124672915181111908317890891659158654221736499522823959933517986673010006749138291836676520080172845444352328"
-        e[7] = big"-0.00514824122639241252178399021479378841872099572255461304439292434131750195489022869965968028106854978547414579491205935930595041763060069987112580994637398395"
-        e[8] = big"0.00296533914055503317169967748114188676589522458557982039693426239853498956125735811263087631479968309978854200615027412311940897061471388689986239742919640848"
-        e[9] = big"-0.0010634368308888065260482548541946175520274736959410047497431569257848032902381738362547705844630238841535652230832162703806430112125115777122361837311714267"
-    else
-        p = num_stages
-        eb = variables(:b, 1:num_stages + 1)
-        @variables y
-        zz = zeros(size(a, 1) + 1)
-        zz2 = zeros(size(a, 1))
-        eA = [zz'
-            zz2 a]
-        ec = [0; c]
-        constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:2*p-3)) do t
-            residual_order_condition(t, RungeKuttaMethod(eA, eb, ec))
-        end 
-        AA, bb, islinear = Symbolics.linear_expansion(Symbolics.substitute.(constraints, (eb[1]=>1/γ,)), eb[2:end])
-        AA = Float64.(map(unwrap, AA))
-        idxs = qr(AA', ColumnNorm()).p[1:num_stages]
-        @assert rank(AA[idxs, :]) == num_stages
-        @assert islinear
-        b_hat = Symbolics.expand.((AA \ -bb))
-        e = [Symbolics.symbolic_to_float(b_hat[i] - b[i]) for i in 1 : num_stages]    
-    end
-    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e, num_stages)
+    # TODO: figure out why all the order conditions are the same
+    A = c_powers'./(1:num_stages)
+    # TODO: figure out why these are the right b
+    b = vcat(-(num_stages)^2, -.5, zeros(num_stages - 2))
+    e = A \ b
+    tab = RadauIIATableau{T1, T2}(T, TI, c2, γ, α, β, e)
 end
+
diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index ab3799c15b..15d7111499 100644
--- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -17,10 +17,10 @@ sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9())
 prob_ode_linear_big = remake(prob_ode_linear, u0 = big.(prob_ode_linear.u0), tspan = big.(prob_ode_linear.tspan))
 prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0), tspan = big.(prob_ode_2Dlinear.tspan))
 
-for i in [3, 5, 7, 9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
+for i in [5, 9, 13], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
     dts = 1 ./ 2 .^ (4.25:-1:0.25)
-    sim21 = test_convergence(dts, prob, AdaptiveRadau(num_stages = i))
-    @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol
+    local sim21 = test_convergence(dts, prob, AdaptiveRadau(min_order = i, max_order = i))
+    @test sim21.𝒪est[:final]≈ i atol=testTol
 end
 
 # test adaptivity
diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_high_order_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_high_order_firk_tests.jl
new file mode 100644
index 0000000000..b78c50f7f1
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_high_order_firk_tests.jl
@@ -0,0 +1,14 @@
+using OrdinaryDiffEqFIRK, DiffEqDevTools, Test, LinearAlgebra
+import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
+
+using GenericSchur
+testTol = 0.5
+
+prob_ode_linear_big = remake(prob_ode_linear, u0 = big.(prob_ode_linear.u0), tspan = big.(prob_ode_linear.tspan))
+prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0), tspan = big.(prob_ode_2Dlinear.tspan))
+
+for i in [17, 21], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
+    dts = 1 ./ 2 .^ (4.25:-1:0.25)
+    sim21 = test_convergence(dts, prob, AdaptiveRadau(min_order = i, max_order = i))
+    @test sim21.𝒪est[:final]≈ i atol=testTol
+end
diff --git a/lib/OrdinaryDiffEqFIRK/test/runtests.jl b/lib/OrdinaryDiffEqFIRK/test/runtests.jl
index 39ede8c3b3..ec335b116b 100644
--- a/lib/OrdinaryDiffEqFIRK/test/runtests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/runtests.jl
@@ -1,3 +1,4 @@
 using SafeTestsets
 
 @time @safetestset "FIRK Tests" include("ode_firk_tests.jl")
+@time @safetestset "High Order FIRK Tests" include("ode_high_order_firk_tests.jl")
diff --git a/lib/OrdinaryDiffEqLinear/Project.toml b/lib/OrdinaryDiffEqLinear/Project.toml
index 107dd521f1..5b8127135a 100644
--- a/lib/OrdinaryDiffEqLinear/Project.toml
+++ b/lib/OrdinaryDiffEqLinear/Project.toml
@@ -18,7 +18,7 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 [compat]
 DiffEqBase = "6.152.2"
 DiffEqDevTools = "2.44.4"
-ExponentialUtilities = "1.26.1"
+ExponentialUtilities = "1.27"
 LinearAlgebra = "<0.0.1, 1"
 OrdinaryDiffEqCore = "1.1"
 OrdinaryDiffEqRosenbrock = "<0.0.1, 1"
diff --git a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl
index 47a1d1bde1..14aaeb0692 100644
--- a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl
+++ b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl
@@ -571,9 +571,9 @@ function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int)
     n = length(u_prototype)
     T = eltype(u_prototype)
     u = zero(u_prototype)                         # stores the current state
-    W = Matrix{T}(undef, n, p + 1)                  # stores the w vectors
-    P = Matrix{T}(undef, n, p + 2)                  # stores output from phiv!
-    Ks = KrylovSubspace{T}(n, maxiter)            # stores output from arnoldi!
+    W = similar(u_prototype, n, p+1)              # stores the w vectors                
+    P = similar(u_prototype, n, p+2)              # stores output from phiv!  
+    Ks = KrylovSubspace{T,T,typeof(similar(u_prototype,size(u_prototype,1),2))}(n, maxiter) # stores output from arnoldi!
     phiv_cache = PhivCache(u_prototype, maxiter, p + 1) # cache used by phiv! (need +1 for error estimation)
     return u, W, P, Ks, phiv_cache
 end
@@ -591,7 +591,7 @@ function alg_cache(alg::LinearExponential, u, rate_prototype, ::Type{uEltypeNoUn
     if alg.krylov == :off
         KsCache = nothing
     elseif alg.krylov == :simple
-        Ks = KrylovSubspace{T}(n, m)
+        Ks = KrylovSubspace{T,T,typeof(similar(u,size(u,1),2))}(n, m)
         expv_cache = ExpvCache{T}(m)
         KsCache = (Ks, expv_cache)
     elseif alg.krylov == :adaptive
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/Project.toml b/lib/OrdinaryDiffEqNonlinearSolve/Project.toml
index 7fdf2ae5c4..27bbe7ff3d 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/Project.toml
+++ b/lib/OrdinaryDiffEqNonlinearSolve/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqNonlinearSolve"
 uuid = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <mayingbo5@gmail.com>"]
-version = "1.2.2"
+version = "1.3.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -35,8 +35,8 @@ ForwardDiff = "0.10.36"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"
 MuladdMacro = "0.2.4"
-NonlinearSolve = "3.14.0"
-OrdinaryDiffEqCore = "1.1"
+NonlinearSolve = "3.14.0, 4"
+OrdinaryDiffEqCore = "1.13"
 OrdinaryDiffEqDifferentiation = "<0.0.1, 1"
 PreallocationTools = "0.4.23"
 Random = "<0.0.1, 1"
@@ -45,7 +45,7 @@ SafeTestsets = "0.1.0"
 SciMLBase = "2.48.1"
 SciMLOperators = "0.3.9"
 SciMLStructures = "1.4.2"
-SimpleNonlinearSolve = "1.12.0"
+SimpleNonlinearSolve = "1.12.0, 2"
 StaticArrays = "1.9.7"
 Test = "<0.0.1, 1"
 julia = "1.10"
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl
index 3453c6fd73..593f27ee77 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl
+++ b/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl
@@ -6,7 +6,8 @@ import SciMLBase
 import SciMLBase: init, solve, solve!, remake
 using SciMLBase: DAEFunction, DEIntegrator, NonlinearFunction, NonlinearProblem,
                  NonlinearLeastSquaresProblem, LinearProblem, ODEProblem, DAEProblem,
-                 update_coefficients!, get_tmp_cache, AbstractSciMLOperator, ReturnCode
+                 update_coefficients!, get_tmp_cache, AbstractSciMLOperator, ReturnCode,
+                 AbstractNonlinearProblem
 import DiffEqBase
 import PreallocationTools
 using SimpleNonlinearSolve: SimpleTrustRegion, SimpleGaussNewton
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl
index 67c3fb177e..b5cffde762 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl
+++ b/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl
@@ -1,5 +1,5 @@
 function default_nlsolve(
-        ::Nothing, isinplace::Val{true}, u, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace::Val{true}, u, ::AbstractNonlinearProblem, autodiff = false)
     FastShortcutNonlinearPolyalg(;
         autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
 end
@@ -8,7 +8,7 @@ function default_nlsolve(
     FastShortcutNLLSPolyalg(; autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
 end
 function default_nlsolve(
-        ::Nothing, isinplace::Val{false}, u, ::NonlinearProblem, autodiff = false)
+        ::Nothing, isinplace::Val{false}, u, ::AbstractNonlinearProblem, autodiff = false)
     FastShortcutNonlinearPolyalg(;
         autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
 end
@@ -17,7 +17,7 @@ function default_nlsolve(
     FastShortcutNLLSPolyalg(; autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
 end
 function default_nlsolve(::Nothing, isinplace::Val{false}, u::StaticArray,
-        ::NonlinearProblem, autodiff = false)
+        ::AbstractNonlinearProblem, autodiff = false)
     SimpleTrustRegion(autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff())
 end
 function default_nlsolve(::Nothing, isinplace::Val{false}, u::StaticArray,
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
index c1c8f67db7..994ecea696 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
+++ b/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
@@ -229,6 +229,10 @@ end
             reltol = reltol)
     end
 
+    if !SciMLBase.successful_retcode(linres.retcode) && linres.retcode != SciMLBase.ReturnCode.Default
+        return convert(eltype(atmp,),Inf)
+    end
+
     cache.linsolve = linres.cache
 
     if DiffEqBase.has_stats(integrator)
diff --git a/lib/OrdinaryDiffEqRosenbrock/Project.toml b/lib/OrdinaryDiffEqRosenbrock/Project.toml
index e837145781..2769329879 100644
--- a/lib/OrdinaryDiffEqRosenbrock/Project.toml
+++ b/lib/OrdinaryDiffEqRosenbrock/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqRosenbrock"
 uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "1.2.0"
+version = "1.3.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl
index d7968cc572..7ed62e5847 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl
@@ -225,7 +225,7 @@ end
 struct RodasTableau{T, T2}
     A::Matrix{T}
     C::Matrix{T}
-    gamma::T
+    gamma::T2
     c::Vector{T2}
     d::Vector{T}
     H::Matrix{T}
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 503484ac33..f6578c3bde 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -58,7 +58,8 @@ import OrdinaryDiffEqCore: trivial_limiter!, CompositeAlgorithm, alg_order,
                            _change_t_via_interpolation!, ODEIntegrator, _ode_interpolant!,
                            current_interpolant, resize_nlsolver!, _ode_interpolant,
                            handle_tstop!, _postamble!, update_uprev!, resize_J_W!,
-                           DAEAlgorithm, get_fsalfirstlast, strip_cache
+                           DAEAlgorithm, get_fsalfirstlast, strip_cache,
+                           Sequential, BaseThreads, PolyesterThreads
 
 export CompositeAlgorithm, ShampineCollocationInit, BrownFullBasicInit, NoInit
 AutoSwitch
diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index 345d90f382..6f934b345f 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -2,7 +2,7 @@
 DDEProblemLibrary = "f42792ee-6ffc-4e2a-ae83-8ee2f22de800"
 DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"	
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
@@ -12,6 +12,7 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 DDEProblemLibrary = "0.1"
 DelayDiffEq = "5.42"
 Measurements = "2.9"
+OrdinaryDiffEq = "6"
 SciMLSensitivity = "7.30"
 StaticArrays = "1"
 StochasticDiffEq = "6.60.1"
diff --git a/test/gpu/linear_exp.jl b/test/gpu/linear_exp.jl
new file mode 100644
index 0000000000..dbbde584cc
--- /dev/null
+++ b/test/gpu/linear_exp.jl
@@ -0,0 +1,34 @@
+using LinearAlgebra
+using SparseArrays
+using CUDA
+using CUDA.CUSPARSE
+using OrdinaryDiffEq
+
+# Linear exponential solvers
+A = MatrixOperator([2.0 -1.0; -1.0 2.0])
+u0 = ones(2)
+
+A_gpu =  MatrixOperator(cu([2.0 -1.0; -1.0 2.0]))
+u0_gpu = cu(ones(2))
+prob_gpu = ODEProblem(A_gpu, u0_gpu, (0.0, 1.0))
+
+sol_analytic = exp(1.0 * Matrix(A)) * u0
+
+@test_broken sol1_gpu = solve(prob_gpu, LinearExponential(krylov = :off))(1.0) |> Vector
+sol2_gpu = solve(prob_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector
+sol3_gpu = solve(prob_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector
+
+@test_broken isapprox(sol1_gpu, sol_analytic, rtol = 1e-6)
+@test isapprox(sol2_gpu, sol_analytic, rtol = 1e-6)
+@test isapprox(sol3_gpu, sol_analytic, rtol = 1e-6)
+
+A2_gpu =  MatrixOperator(cu(sparse([2.0 -1.0; -1.0 2.0])))
+prob2_gpu = ODEProblem(A2_gpu, u0_gpu, (0.0, 1.0))
+
+@test_broken sol2_1_gpu = solve(prob2_gpu, LinearExponential(krylov = :off))(1.0) |> Vector
+sol2_2_gpu = solve(prob2_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector
+sol2_3_gpu = solve(prob2_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector
+
+@test_broken isapprox(sol2_1_gpu, sol_analytic, rtol = 1e-6)
+@test isapprox(sol2_2_gpu, sol_analytic, rtol = 1e-6)
+@test isapprox(sol2_3_gpu, sol_analytic, rtol = 1e-6)
\ No newline at end of file
diff --git a/test/interface/checkinit_tests.jl b/test/interface/checkinit_tests.jl
index 8320449151..e71b412e2c 100644
--- a/test/interface/checkinit_tests.jl
+++ b/test/interface/checkinit_tests.jl
@@ -24,9 +24,9 @@ roberf_oop = ODEFunction{false}(rober, mass_matrix = M)
 prob_mm = ODEProblem(roberf, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4))
 prob_mm_oop = ODEProblem(roberf_oop, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4))
 
-@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(
+@test_throws SciMLBase.CheckInitFailureError solve(
     prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit())
-@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(
+@test_throws SciMLBase.CheckInitFailureError solve(
     prob_mm_oop, Rodas5P(), reltol = 1e-8, abstol = 1e-8,
     initializealg = SciMLBase.CheckInit())
 
@@ -49,7 +49,7 @@ tspan = (0.0, 100000.0)
 differential_vars = [true, true, false]
 prob = DAEProblem(f, du₀, u₀, tspan, differential_vars = differential_vars)
 prob_oop = DAEProblem(f_oop, du₀, u₀, tspan, differential_vars = differential_vars)
-@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(
+@test_throws SciMLBase.CheckInitFailureError solve(
     prob, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit())
-@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(
+@test_throws SciMLBase.CheckInitFailureError solve(
     prob_oop, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit())
diff --git a/test/interface/ode_tstops_tests.jl b/test/interface/ode_tstops_tests.jl
index ce85f0e859..a911ecfe8a 100644
--- a/test/interface/ode_tstops_tests.jl
+++ b/test/interface/ode_tstops_tests.jl
@@ -76,3 +76,15 @@ end
     prob = ODEProblem(ff, [0.0], (0.0f0, 1.0f0))
     sol = solve(prob, Tsit5(), tstops = [tval], callback = cb)
 end
+
+@testset "Late binding tstops" begin
+    function rhs(u, p, t)
+        u * p + t
+    end
+    prob = ODEProblem(rhs, 1.0, (0.0, 1.0), 0.1; tstops = (p, tspan) -> tspan[1]:p:tspan[2])
+    sol = solve(prob, Tsit5())
+    @test 0.0:0.1:1.0 ⊆ sol.t
+    prob2 = remake(prob; p = 0.07)
+    sol2 = solve(prob2, Tsit5())
+    @test 0.0:0.07:1.0 ⊆ sol2.t
+end
diff --git a/test/interface/precision_mixing.jl b/test/interface/precision_mixing.jl
new file mode 100644
index 0000000000..497944d65d
--- /dev/null
+++ b/test/interface/precision_mixing.jl
@@ -0,0 +1,22 @@
+using OrdinaryDiffEq
+function ODE(du, u, t, R, K)
+    du .= u
+end
+params = BigFloat[1. 0.91758707304098; 1.48439909482661 1.]
+u0 = BigFloat[0.1, 0.1]
+tspan = (1.0, 31.0)
+R = BigFloat[0.443280390004304303, 1.172917082211452]
+K = BigFloat[13.470600276901400604, 12.52980757005]
+ODE_ = (du, u, params, t) -> ODE(du, u, t, R, K)
+odeProblem = ODEProblem(ODE_, u0, tspan, params)
+for alg in [AutoVern8(Rodas5(), nonstifftol = 11 / 10)
+            FBDF()
+            QNDF()
+            Tsit5()
+            Rodas5P()
+            TRBDF2()
+            KenCarp4()
+            RadauIIA5()
+    ]
+    Solution = solve(odeProblem, alg, saveat = 1, abstol = 1.e-12, reltol = 1.e-6)
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 821c48615d..55b58486ab 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -85,6 +85,7 @@ end
     if !is_APPVEYOR && (GROUP == "All" || GROUP == "InterfaceIV" || GROUP == "Interface")
         @time @safetestset "Autodiff Error Tests" include("interface/autodiff_error_tests.jl")
         @time @safetestset "Ambiguity Tests" include("interface/ambiguity_tests.jl")
+        @time @safetestset "Precision Mixing Tests" include("interface/precision_mixing.jl")
         @time @safetestset "Sized Matrix Tests" include("interface/sized_matrix_tests.jl")
         @time @safetestset "Second Order with First Order Solver Tests" include("interface/second_order_with_first_order_solvers.jl")
     end
@@ -171,6 +172,7 @@ end
         end
         @time @safetestset "Autoswitch GPU" include("gpu/autoswitch.jl")
         @time @safetestset "Linear LSRK GPU" include("gpu/linear_lsrk.jl")
+        @time @safetestset "Linear Exponential GPU" include("gpu/linear_exp.jl")
         @time @safetestset "Reaction-Diffusion Stiff Solver GPU" include("gpu/reaction_diffusion_stiff.jl")
         @time @safetestset "Scalar indexing bug bypass" include("gpu/hermite_test.jl")
     end