diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 78cb0d9372..54ad492b93 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -34,6 +34,7 @@ jobs:
           - OrdinaryDiffEqExponentialRK
           - OrdinaryDiffEqExtrapolation
           - OrdinaryDiffEqFIRK
+          - OrdinaryDiffEqFIRKGenerator
           - OrdinaryDiffEqFeagin
           - OrdinaryDiffEqFunctionMap
           - OrdinaryDiffEqHighOrderRK
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/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..ea0eed4466 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.11.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.60"
 SciMLOperators = "0.3"
 SciMLStructures = "1"
 SimpleUnPack = "1"
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/initialize_dae.jl b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl
index e88deeef1a..ebefa3d91c 100644
--- a/lib/OrdinaryDiffEqCore/src/initialize_dae.jl
+++ b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl
@@ -169,7 +169,7 @@ function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
     end
 
     alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD)
-    nlsol = solve(initializeprob, alg)
+    nlsol = solve(initializeprob, alg, abstol = integrator.opts.abstol, reltol = integrator.opts.reltol)
     if isinplace === Val{true}()
         integrator.u .= prob.f.initializeprobmap(nlsol)
     elseif isinplace === Val{false}()
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/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..0525c0059e 100644
--- a/lib/OrdinaryDiffEqFIRK/Project.toml
+++ b/lib/OrdinaryDiffEqFIRK/Project.toml
@@ -1,34 +1,27 @@
 name = "OrdinaryDiffEqFIRK"
 uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "1.2.0"
+version = "1.4.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 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"
 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"
 FastPower = "1"
-GenericLinearAlgebra = "0.3.13"
-GenericSchur = "0.5.4"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"
 MuladdMacro = "0.2.4"
@@ -36,14 +29,11 @@ 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"
 SciMLOperators = "0.3.9"
-Symbolics = "6.15.3"
 Test = "<0.0.1, 1"
 julia = "1.10"
 
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/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index f0684fb259..99374e959f 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -528,125 +528,6 @@ function BigRadauIIA13Tableau(T1, T2)
     c, γ, α, β, e)
 end
 
-using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
-using Symbolics: variables, variable, unwrap
-
-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)
-    end
-    c = real(roots(radau_p))
-    c[num_stages] = 1
-    c_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
-    for i in 1 : num_stages
-        for j in 1 : num_stages
-            c_powers[i,j] = c[i]^(j - 1)
-        end
-    end
-    inverse_c_powers = inv(c_powers)
-    c_q = Matrix{BigFloat}(undef, num_stages, num_stages)
-    for i in 1 : num_stages
-        for j in 1 : num_stages
-            c_q[i,j] = c[i]^(j) / j
-        end
-    end
-    a = c_q * inverse_c_powers
-    a_inverse = inv(a)
-    b = Vector{BigFloat}(undef, 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
-    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
-    end
-    TI = inv(T)
-
-    if (num_stages == 9)
-        e = Vector{BigFloat}(undef, 9)
-        e[1] = big"-89.8315397040376845865027298766511166861131537901479318008187013574099993398844876573472315778350373191126204142357525815115482293843777624541394691345885716"
-        e[2] = big"11.4742766094687721590222610299234578063148408248968597722844661019124491691448775794163842022854672278004372474682761156236829237591471118886342174262239472"
-        e[3] = big"-3.81419058476042873698615187248837320040477891376179026064712181641592908409919668221598902628694008903410444392769866137859041139561191341971835412426311966"
-        e[4] = big"1.81155300867853110911564243387531599775142729190474576183505286509346678884073482369609308584446518479366940471952219053256362416491879701351428578466580598"
-        e[5] = big"-1.03663781378817415276482837566889343026914084945266083480559060702535168750966084568642219911350874500410428043808038021858812311835772945467924877281164517"
-        e[6] = big"0.660865688193716483757690045578935452512421753840843511309717716369201467579470723336314286637650332622546110594223451602017981477424498704954672224534648119"
-        e[7] = big"-0.444189256280526730087023435911479370800996444567516110958885112499737452734669537494435549195615660656770091500773942469075264796140815048410568498349675229"
-        e[8] = big"0.290973163636905565556251162453264542120491238398561072912173321087011249774042707406397888774630179702057578431394918930648610404108923880955576205699885598"
-        e[9] = big"-0.111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111222795"
-    elseif (num_stages == 11)
-        e = Vector{BigFloat}(undef, 11)
-        e[1] = big"-134.152626015465044063378550835075318643291579891352838474367124350171545245813244797505763447327562609902792066283575334085390478517120485782603677022267543"
-        e[2] = big"17.0660253399060146849212356299749772423073416838121578997449942694355150369717420038613850964748566731121793290881077515821557030349184664685171028112845693"
-        e[3] = big"-5.63464089555106294823267450977601185069165875295372865523759287935369597689662768988715406731927279137711764532851201746616033935275093116699140897901326857"
-        e[4] = big"2.65398285960564394428637524662555134392389271086844331137910389226095922845489762567700560496915255196379049844894623384211693438658842276927416827629120392"
-        e[5] = big"-1.50753272514563441873424939425410006034401178578882643601844794171149654717227697249290904230103304153661631200445957060050895700394738491883951084826421405"
-        e[6] = big"0.960260572218344245935269463733859188992760928707230734981795807797858324380878500135029848170473080912207529262984056182004711806457345405466997261506487216"
-        e[7] = big"-0.658533932484491373507110339620843007350146695468297825313721271556868110859353953892288534787571420691760379406525738632649863532050280264983313133523641674"
-        e[8] = big"0.47189364490739958527881800092758816959227958959727295348380187162217987951960275929676019062173412149363239153353720640122975284789262792027244826613784432"
-        e[9] = big"-0.34181016557091711933253384050957887606039737751222218385118573305954222606860932803075900338195356026497059819558648780544900376040113065955083806288937526"
-        e[10] = big"0.233890408488838371854329668882967402012428680999899584289285425645726546573900943747784263972086087200538161975992991491742449181322441138528940521648041699"
-        e[11] = big"-0.0909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909093788951"
-    elseif (num_stages == 13)
-        e = Vector{BigFloat}(undef, 13)
-        e[1] = big"-187.337806666035250696387113105488477375830948862159770885826492736743460038872636916422100706359786154665214547894636085276885830138994748219148357620227002"
-        e[2] = big"23.775705048946302520021716862887025159493544949407763131913924588605891085865877529749667170060976683489861224477421212170329019074926368036881685518012728"
-        e[3] = big"-7.81823724708755833325842676798052630403951326380926053607036280237871312516353176794790424805918285990907426633641064901501063343970205708057561515795364672"
-        e[4] = big"3.66289388251066047904501665386587373682645522696191680651425553890800106379174431775463608296821504040006089759980653462003322200870566661322334735061646223"
-        e[5] = big"-2.06847094952801462392548700163367193433237251061765813625197254100990426184032443671875204952150187523419743001493620194301209589692419776688692360679336566"
-        e[6] = big"1.31105635982993157063104433803023633257356281733787535204132865785504258558244947718491624714070193102812968996631302993877989767202703509685785407541965509"
-        e[7] = big"-0.897988270828178667954874573865888835427640297795141000639881363403080887358272161865529150995401606679722232843051402663087372891040498351714982629218397165"
-        e[8] = big"0.648958340079591709325028357505725843500310779765000237611355105578356380892509437805732950287939403489669590070670546599339082534053791877148407548785389408"
-        e[9] = big"-0.485906120880156534303797908584178831869407602334908394589833216071089678420073112977712585616439120156658051446412515753614726507868506301824972455936531663"
-        e[10] = big"0.370151313405058266144090771980402238126294149688261261935258556082315591034906662511634673912342573394958760869036835172495369190026354174118335052418701339"
-        e[11] = big"-0.27934271062931554435643589252670994638477019847143394253283050767117135003630906657393675748475838251860910095199485920686192935009874559019443503474805827"
-        e[12] = big"0.195910097140006778096161342733266840441407888950433028972173797170889557600583114422425296743817444283872389581116632280572920821812614435192580036549169031"
-        e[13] = big"-0.0769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769254590189"
-    else
-        e_sym = variables(:e, 1:num_stages)
-        constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:num_stages)) do t
-            residual_order_condition(t, RungeKuttaMethod(a, e_sym, c))
-        end 
-        AA, bb, islinear = Symbolics.linear_expansion(constraints, e_sym[1:end])
-        AA = BigFloat.(map(unwrap, AA))
-        bb = BigFloat.(map(unwrap, bb))
-        A = vcat([zeros(num_stages -1); 1]', AA)
-        b_2 = vcat(-1/big(num_stages), -(num_stages)^2, -1, zeros(size(A, 1) - 3))
-        e = A \ b_2
-    end
-    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e)
+function adaptiveRadauTableau(T1, T2, num_stages)
+    error("num_stages choice $num_stages out of the pre-generated tableau range. For the fully adaptive Radau, please load the extension via `using OrdinaryDiffEqFIRKGenerator`")
 end
diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index ad479d172d..72e0be40f4 100644
--- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -17,7 +17,7 @@ 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 [3, 5, 7], 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_stages = i, max_stages = i))
     @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol
diff --git a/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md b/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md
new file mode 100644
index 0000000000..4a7df96ac5
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRKGenerator/LICENSE.md
@@ -0,0 +1,24 @@
+The OrdinaryDiffEq.jl package is licensed under the MIT "Expat" License:
+
+> Copyright (c) 2016-2020: ChrisRackauckas, Yingbo Ma, Julia Computing Inc, and
+> other contributors:
+> 
+> https://github.com/SciML/OrdinaryDiffEq.jl/graphs/contributors
+> 
+> Permission is hereby granted, free of charge, to any person obtaining a copy
+> of this software and associated documentation files (the "Software"), to deal
+> in the Software without restriction, including without limitation the rights
+> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+> copies of the Software, and to permit persons to whom the Software is
+> furnished to do so, subject to the following conditions:
+> 
+> The above copyright notice and this permission notice shall be included in all
+> copies or substantial portions of the Software.
+> 
+> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+> SOFTWARE.
diff --git a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml
new file mode 100644
index 0000000000..3e4156ddea
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml
@@ -0,0 +1,35 @@
+name = "OrdinaryDiffEqFIRKGenerator"
+uuid = "677d4f02-548a-44fa-8eaf-26579094acaf"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.0.0"
+
+[deps]
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrdinaryDiffEqFIRK = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+
+[compat]
+DiffEqDevTools = "2.44.4"
+GenericLinearAlgebra = "0.3.13"
+GenericSchur = "0.5.4"
+LinearAlgebra = "<0.0.1, 1"
+ODEProblemLibrary = "0.1.8"
+OrdinaryDiffEqFIRK = "1"
+Polynomials = "4.0.11"
+RootedTrees = "2.23.1"
+Symbolics = "6.15.3"
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+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"]
diff --git a/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl b/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl
new file mode 100644
index 0000000000..f205bf47ce
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRKGenerator/src/OrdinaryDiffEqFIRKGenerator.jl
@@ -0,0 +1,127 @@
+module OrdinaryDiffEqFIRKGenerator
+
+using OrdinaryDiffEqFIRK
+using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
+using Symbolics: variables, variable, unwrap
+
+function OrdinaryDiffEqFIRK.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)
+    end
+    c = real(roots(radau_p))
+    c[num_stages] = 1
+    c_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
+    for i in 1 : num_stages
+        for j in 1 : num_stages
+            c_powers[i,j] = c[i]^(j - 1)
+        end
+    end
+    inverse_c_powers = inv(c_powers)
+    c_q = Matrix{BigFloat}(undef, num_stages, num_stages)
+    for i in 1 : num_stages
+        for j in 1 : num_stages
+            c_q[i,j] = c[i]^(j) / j
+        end
+    end
+    a = c_q * inverse_c_powers
+    a_inverse = inv(a)
+    b = Vector{BigFloat}(undef, 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
+    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
+    end
+    TI = inv(T)
+
+    if (num_stages == 9)
+        e = Vector{BigFloat}(undef, 9)
+        e[1] = big"-89.8315397040376845865027298766511166861131537901479318008187013574099993398844876573472315778350373191126204142357525815115482293843777624541394691345885716"
+        e[2] = big"11.4742766094687721590222610299234578063148408248968597722844661019124491691448775794163842022854672278004372474682761156236829237591471118886342174262239472"
+        e[3] = big"-3.81419058476042873698615187248837320040477891376179026064712181641592908409919668221598902628694008903410444392769866137859041139561191341971835412426311966"
+        e[4] = big"1.81155300867853110911564243387531599775142729190474576183505286509346678884073482369609308584446518479366940471952219053256362416491879701351428578466580598"
+        e[5] = big"-1.03663781378817415276482837566889343026914084945266083480559060702535168750966084568642219911350874500410428043808038021858812311835772945467924877281164517"
+        e[6] = big"0.660865688193716483757690045578935452512421753840843511309717716369201467579470723336314286637650332622546110594223451602017981477424498704954672224534648119"
+        e[7] = big"-0.444189256280526730087023435911479370800996444567516110958885112499737452734669537494435549195615660656770091500773942469075264796140815048410568498349675229"
+        e[8] = big"0.290973163636905565556251162453264542120491238398561072912173321087011249774042707406397888774630179702057578431394918930648610404108923880955576205699885598"
+        e[9] = big"-0.111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111222795"
+    elseif (num_stages == 11)
+        e = Vector{BigFloat}(undef, 11)
+        e[1] = big"-134.152626015465044063378550835075318643291579891352838474367124350171545245813244797505763447327562609902792066283575334085390478517120485782603677022267543"
+        e[2] = big"17.0660253399060146849212356299749772423073416838121578997449942694355150369717420038613850964748566731121793290881077515821557030349184664685171028112845693"
+        e[3] = big"-5.63464089555106294823267450977601185069165875295372865523759287935369597689662768988715406731927279137711764532851201746616033935275093116699140897901326857"
+        e[4] = big"2.65398285960564394428637524662555134392389271086844331137910389226095922845489762567700560496915255196379049844894623384211693438658842276927416827629120392"
+        e[5] = big"-1.50753272514563441873424939425410006034401178578882643601844794171149654717227697249290904230103304153661631200445957060050895700394738491883951084826421405"
+        e[6] = big"0.960260572218344245935269463733859188992760928707230734981795807797858324380878500135029848170473080912207529262984056182004711806457345405466997261506487216"
+        e[7] = big"-0.658533932484491373507110339620843007350146695468297825313721271556868110859353953892288534787571420691760379406525738632649863532050280264983313133523641674"
+        e[8] = big"0.47189364490739958527881800092758816959227958959727295348380187162217987951960275929676019062173412149363239153353720640122975284789262792027244826613784432"
+        e[9] = big"-0.34181016557091711933253384050957887606039737751222218385118573305954222606860932803075900338195356026497059819558648780544900376040113065955083806288937526"
+        e[10] = big"0.233890408488838371854329668882967402012428680999899584289285425645726546573900943747784263972086087200538161975992991491742449181322441138528940521648041699"
+        e[11] = big"-0.0909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909093788951"
+    elseif (num_stages == 13)
+        e = Vector{BigFloat}(undef, 13)
+        e[1] = big"-187.337806666035250696387113105488477375830948862159770885826492736743460038872636916422100706359786154665214547894636085276885830138994748219148357620227002"
+        e[2] = big"23.775705048946302520021716862887025159493544949407763131913924588605891085865877529749667170060976683489861224477421212170329019074926368036881685518012728"
+        e[3] = big"-7.81823724708755833325842676798052630403951326380926053607036280237871312516353176794790424805918285990907426633641064901501063343970205708057561515795364672"
+        e[4] = big"3.66289388251066047904501665386587373682645522696191680651425553890800106379174431775463608296821504040006089759980653462003322200870566661322334735061646223"
+        e[5] = big"-2.06847094952801462392548700163367193433237251061765813625197254100990426184032443671875204952150187523419743001493620194301209589692419776688692360679336566"
+        e[6] = big"1.31105635982993157063104433803023633257356281733787535204132865785504258558244947718491624714070193102812968996631302993877989767202703509685785407541965509"
+        e[7] = big"-0.897988270828178667954874573865888835427640297795141000639881363403080887358272161865529150995401606679722232843051402663087372891040498351714982629218397165"
+        e[8] = big"0.648958340079591709325028357505725843500310779765000237611355105578356380892509437805732950287939403489669590070670546599339082534053791877148407548785389408"
+        e[9] = big"-0.485906120880156534303797908584178831869407602334908394589833216071089678420073112977712585616439120156658051446412515753614726507868506301824972455936531663"
+        e[10] = big"0.370151313405058266144090771980402238126294149688261261935258556082315591034906662511634673912342573394958760869036835172495369190026354174118335052418701339"
+        e[11] = big"-0.27934271062931554435643589252670994638477019847143394253283050767117135003630906657393675748475838251860910095199485920686192935009874559019443503474805827"
+        e[12] = big"0.195910097140006778096161342733266840441407888950433028972173797170889557600583114422425296743817444283872389581116632280572920821812614435192580036549169031"
+        e[13] = big"-0.0769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769254590189"
+    else
+        e_sym = variables(:e, 1:num_stages)
+        constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:num_stages)) do t
+            residual_order_condition(t, RungeKuttaMethod(a, e_sym, c))
+        end 
+        AA, bb, islinear = Symbolics.linear_expansion(constraints, e_sym[1:end])
+        AA = BigFloat.(map(unwrap, AA))
+        bb = BigFloat.(map(unwrap, bb))
+        A = vcat([zeros(num_stages -1); 1]', AA)
+        b_2 = vcat(-1/big(num_stages), -(num_stages)^2, -1, zeros(size(A, 1) - 3))
+        e = A \ b_2
+    end
+    OrdinaryDiffEqFIRK.RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e)
+end
+
+end
diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl
new file mode 100644
index 0000000000..d3c19ebbb5
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl
@@ -0,0 +1,13 @@
+using OrdinaryDiffEqFIRK, OrdinaryDiffEqFIRKGenerator, DiffEqDevTools, Test, LinearAlgebra
+import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
+
+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 [9], 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_stages = i, max_stages = i))
+    @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol
+end
diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl
new file mode 100644
index 0000000000..108f9267b9
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRKGenerator/test/runtests.jl
@@ -0,0 +1,3 @@
+using SafeTestsets
+
+@time @safetestset "Generated FIRK Tests" include("ode_firk_tests.jl")
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/Project.toml b/lib/OrdinaryDiffEqNonlinearSolve/Project.toml
index 7fdf2ae5c4..68f3db0d8e 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.2.4"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -35,7 +35,7 @@ ForwardDiff = "0.10.36"
 LinearAlgebra = "<0.0.1, 1"
 LinearSolve = "2.32.0"
 MuladdMacro = "0.2.4"
-NonlinearSolve = "3.14.0"
+NonlinearSolve = "3.14.0, 4"
 OrdinaryDiffEqCore = "1.1"
 OrdinaryDiffEqDifferentiation = "<0.0.1, 1"
 PreallocationTools = "0.4.23"
@@ -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/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/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..3928ddad86 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