diff --git a/lib/OrdinaryDiffEqExtrapolation/test/runtests.jl b/lib/OrdinaryDiffEqExtrapolation/test/runtests.jl
index 738c8013cc..4aee75dc86 100644
--- a/lib/OrdinaryDiffEqExtrapolation/test/runtests.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/test/runtests.jl
@@ -1,3 +1,3 @@
 using SafeTestsets
 
-@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")
+@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")
\ No newline at end of file
diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml
new file mode 100644
index 0000000000..e87c1dd521
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/Project.toml
@@ -0,0 +1,24 @@
+name = "OrdinaryDiffEqFIRK"
+uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.0.0"
+
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[compat]
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]
\ No newline at end of file
diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
new file mode 100644
index 0000000000..d49edd4a1b
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
@@ -0,0 +1,35 @@
+module OrdinaryDiffEqFIRK
+
+import OrdinaryDiffEq: alg_order, calculate_residuals!,
+                       initialize!, perform_step!, @unpack, unwrap_alg,
+                       calculate_residuals,
+                       OrdinaryDiffEqAlgorithm, OrdinaryDiffEqNewtonAdaptiveAlgorithm,
+                       OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
+                       OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
+                       alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
+                       constvalue, _unwrap_val, du_alias_or_new,
+                       explicit_rk_docstring, trivial_limiter!,
+                       _ode_interpolant!, _ode_addsteps!, AbstractController,
+                       NLStatus, qmax_default, alg_adaptive_order, DEFAULT_PRECS,
+                       UJacobianWrapper, build_J_W, build_jac_config, UDerivativeWrapper,
+                       Convergence, calc_J!, dolinsolve, FastConvergence, calc_J,
+                       stepsize_controller!, islinearfunction, step_accept_controller!, step_reject_controller!,
+                       PredictiveController, alg_can_repeat_jac, NewtonAlgorithm, fac_default_gamma,
+                       get_new_W_γdt_cutoff, get_current_adaptive_order, VerySlowConvergence, 
+                       Divergence, isfirk
+using MuladdMacro, DiffEqBase, RecursiveArrayTools
+using SciMLOperators: AbstractSciMLOperator
+using LinearAlgebra: I, UniformScaling, mul!, lu
+import LinearSolve
+import FastBroadcast: @..
+include("algorithms.jl")
+include("alg_utils.jl")
+include("controllers.jl")
+include("firk_caches.jl")
+include("firk_tableaus.jl")
+include("firk_perform_step.jl")
+include("integrator_interface.jl")
+
+export RadauIIA3, RadauIIA5, RadauIIA7
+
+end
diff --git a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
new file mode 100644
index 0000000000..ae8acba221
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
@@ -0,0 +1,13 @@
+qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA7}) = 8
+
+alg_order(alg::RadauIIA3) = 3
+alg_order(alg::RadauIIA5) = 5
+alg_order(alg::RadauIIA7) = 7
+
+isfirk(alg::RadauIIA3) = true
+isfirk(alg::RadauIIA5) = true
+isfirk(alg::RadauIIA7) = true
+
+alg_adaptive_order(alg::RadauIIA3) = 1
+alg_adaptive_order(alg::RadauIIA5) = 3
+alg_adaptive_order(alg::RadauIIA7) = 5
diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
new file mode 100644
index 0000000000..05953b9d58
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -0,0 +1,155 @@
+# FIRK Methods
+
+"""
+@article{hairer1999stiff,
+title={Stiff differential equations solved by Radau methods},
+author={Hairer, Ernst and Wanner, Gerhard},
+journal={Journal of Computational and Applied Mathematics},
+volume={111},
+number={1-2},
+pages={93--111},
+year={1999},
+publisher={Elsevier}
+}
+
+RadauIIA3: Fully-Implicit Runge-Kutta Method
+An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
+"""
+struct RadauIIA3{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    precs::P
+    extrapolant::Symbol
+    κ::Tol
+    maxiters::Int
+    fast_convergence_cutoff::C1
+    new_W_γdt_cutoff::C2
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function RadauIIA3(; chunk_size = Val{0}(), autodiff = Val{true}(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS,
+        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+        new_W_γdt_cutoff = 1 // 5,
+        controller = :Predictive, κ = nothing, maxiters = 10,
+        step_limiter! = trivial_limiter!)
+    RadauIIA3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(fast_convergence_cutoff),
+        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+        precs,
+        extrapolant,
+        κ,
+        maxiters,
+        fast_convergence_cutoff,
+        new_W_γdt_cutoff,
+        controller,
+        step_limiter!)
+end
+
+"""
+@article{hairer1999stiff,
+title={Stiff differential equations solved by Radau methods},
+author={Hairer, Ernst and Wanner, Gerhard},
+journal={Journal of Computational and Applied Mathematics},
+volume={111},
+number={1-2},
+pages={93--111},
+year={1999},
+publisher={Elsevier}
+}
+
+RadauIIA5: Fully-Implicit Runge-Kutta Method
+An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
+"""
+struct RadauIIA5{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    precs::P
+    smooth_est::Bool
+    extrapolant::Symbol
+    κ::Tol
+    maxiters::Int
+    fast_convergence_cutoff::C1
+    new_W_γdt_cutoff::C2
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function RadauIIA5(; chunk_size = Val{0}(), autodiff = Val{true}(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS,
+        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+        new_W_γdt_cutoff = 1 // 5,
+        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
+        step_limiter! = trivial_limiter!)
+    RadauIIA5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(fast_convergence_cutoff),
+        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+        precs,
+        smooth_est,
+        extrapolant,
+        κ,
+        maxiters,
+        fast_convergence_cutoff,
+        new_W_γdt_cutoff,
+        controller,
+        step_limiter!)
+end
+
+"""
+@article{hairer1999stiff,
+title={Stiff differential equations solved by Radau methods},
+author={Hairer, Ernst and Wanner, Gerhard},
+journal={Journal of Computational and Applied Mathematics},
+volume={111},
+number={1-2},
+pages={93--111},
+year={1999},
+publisher={Elsevier}
+}
+
+RadauIIA7: Fully-Implicit Runge-Kutta Method
+An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
+"""
+struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    precs::P
+    smooth_est::Bool
+    extrapolant::Symbol
+    κ::Tol
+    maxiters::Int
+    fast_convergence_cutoff::C1
+    new_W_γdt_cutoff::C2
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function RadauIIA7(; chunk_size = Val{0}(), autodiff = Val{true}(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS,
+        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+        new_W_γdt_cutoff = 1 // 5,
+        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
+        step_limiter! = trivial_limiter!)
+    RadauIIA7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(fast_convergence_cutoff),
+        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+        precs,
+        smooth_est,
+        extrapolant,
+        κ,
+        maxiters,
+        fast_convergence_cutoff,
+        new_W_γdt_cutoff,
+        controller,
+        step_limiter!)
+end
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
new file mode 100644
index 0000000000..df6cf153ab
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -0,0 +1,51 @@
+@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 = DiffEqBase.fastpow(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) *
+               DiffEqBase.fastpow((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
diff --git a/src/caches/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
similarity index 99%
rename from src/caches/firk_caches.jl
rename to lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index fb1489270b..5a560d2ae4 100644
--- a/src/caches/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -200,7 +200,6 @@ mutable struct RadauIIA5Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     status::NLStatus
     step_limiter!::StepLimiter
 end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA5Cache 1
 
 function alg_cache(alg::RadauIIA5, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits},
@@ -370,7 +369,6 @@ mutable struct RadauIIA7Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     status::NLStatus
     step_limiter!::StepLimiter
 end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA7Cache 1
 
 function alg_cache(alg::RadauIIA7, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits},
diff --git a/src/perform_step/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
similarity index 100%
rename from src/perform_step/firk_perform_step.jl
rename to lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
diff --git a/src/tableaus/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
similarity index 100%
rename from src/tableaus/firk_tableaus.jl
rename to lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
diff --git a/lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl b/lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl
new file mode 100644
index 0000000000..5266edfdc8
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl
@@ -0,0 +1,4 @@
+@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA7},
+        cache::OrdinaryDiffEqMutableCache)
+    (cache.tmp, cache.atmp)
+end
diff --git a/test/algconvergence/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
similarity index 99%
rename from test/algconvergence/ode_firk_tests.jl
rename to lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index 3672e30fa9..a7003ce88a 100644
--- a/test/algconvergence/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -57,4 +57,4 @@ for iip in (true, false)
         @test sol.stats.njacs < sol.stats.nw # W reuse
     end
     @test length(sol) < 5000 # the error estimate is not very good
-end
+end
\ No newline at end of file
diff --git a/lib/OrdinaryDiffEqFIRK/test/runtests.jl b/lib/OrdinaryDiffEqFIRK/test/runtests.jl
new file mode 100644
index 0000000000..b252cdf8f2
--- /dev/null
+++ b/lib/OrdinaryDiffEqFIRK/test/runtests.jl
@@ -0,0 +1,3 @@
+using SafeTestsets
+
+@time @safetestset "FIRK Tests" include("ode_firk_tests.jl")
\ No newline at end of file
diff --git a/lib/OrdinaryDiffEqFeagin/test/runtests.jl b/lib/OrdinaryDiffEqFeagin/test/runtests.jl
index 8b13789179..e24bef2c88 100644
--- a/lib/OrdinaryDiffEqFeagin/test/runtests.jl
+++ b/lib/OrdinaryDiffEqFeagin/test/runtests.jl
@@ -1 +1,4 @@
 
+using SafeTestsets
+
+@time @safetestset "Feagin Tests" include("ode_feagin_tests.jl")
\ No newline at end of file
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 6214e87808..362e21ff58 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -148,7 +148,6 @@ include("generic_rosenbrock.jl")
 include("caches/basic_caches.jl")
 include("caches/low_order_rk_caches.jl")
 include("caches/high_order_rk_caches.jl")
-include("caches/firk_caches.jl")
 include("caches/linear_caches.jl")
 include("caches/linear_nonlinear_caches.jl")
 include("caches/rosenbrock_caches.jl")
@@ -162,7 +161,6 @@ include("caches/qprk_caches.jl")
 include("tableaus/low_order_rk_tableaus.jl")
 include("tableaus/high_order_rk_tableaus.jl")
 include("tableaus/rosenbrock_tableaus.jl")
-include("tableaus/firk_tableaus.jl")
 include("tableaus/qprk_tableaus.jl")
 
 include("integrators/type.jl")
@@ -181,7 +179,6 @@ include("perform_step/exponential_rk_perform_step.jl")
 include("perform_step/explicit_rk_perform_step.jl")
 include("perform_step/low_order_rk_perform_step.jl")
 include("perform_step/high_order_rk_perform_step.jl")
-include("perform_step/firk_perform_step.jl")
 include("perform_step/rosenbrock_perform_step.jl")
 include("perform_step/composite_perform_step.jl")
 include("perform_step/adams_bashforth_moulton_perform_step.jl")
@@ -284,11 +281,16 @@ include("../lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl")
 using ..OrdinaryDiffEqDefault
 export DefaultODEAlgorithm
 
+include("../lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl")
+using ..OrdinaryDiffEqFIRK
+export RadauIIA3, RadauIIA5, RadauIIA7
+
 using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!,
                            calc_Lagrange_interp!,
                            calc_finite_difference_weights, estimate_terk,
                            calc_Lagrange_interp,
                            bdf_step_reject_controller!
+                           
 include("nlsolve/newton.jl")
 include("perform_step/dae_perform_step.jl")
 
@@ -428,8 +430,6 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
        FRK65, PFRK87,
        RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4
 
-export RadauIIA3, RadauIIA5, RadauIIA7
-
 export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
        MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4,
        MagnusAdapt4, RKMK2, RKMK4, LieRK4, CG2, CG3, CG4a
diff --git a/src/alg_utils.jl b/src/alg_utils.jl
index 77bf1df0f4..4120c7652b 100644
--- a/src/alg_utils.jl
+++ b/src/alg_utils.jl
@@ -50,6 +50,8 @@ isfsal(alg::PSRK3p5q4) = false
 isfsal(alg::PSRK3p6q5) = false
 isfsal(alg::PSRK4p7q6) = false
 
+isfirk(alg) = false
+
 get_current_isfsal(alg, cache) = isfsal(alg)
 
 # evaluates f(t[i])
@@ -175,7 +177,6 @@ qmin_default(alg::DP8) = 1 // 3
 qmax_default(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = 10
 qmax_default(alg::CompositeAlgorithm) = minimum(qmax_default.(alg.algs))
 qmax_default(alg::DP8) = 6
-qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA7}) = 8
 
 function has_chunksize(alg::OrdinaryDiffEqAlgorithm)
     return alg isa Union{OrdinaryDiffEqExponentialAlgorithm,
@@ -428,9 +429,6 @@ alg_order(alg::Tsit5) = 5
 alg_order(alg::DP8) = 8
 alg_order(alg::TanYam7) = 7
 alg_order(alg::TsitPap8) = 8
-alg_order(alg::RadauIIA3) = 3
-alg_order(alg::RadauIIA5) = 5
-alg_order(alg::RadauIIA7) = 7
 alg_order(alg::RKMK2) = 2
 alg_order(alg::RKMK4) = 4
 alg_order(alg::LieRK4) = 4
@@ -528,9 +526,6 @@ alg_adaptive_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = alg_orde
 alg_adaptive_order(alg::Rosenbrock23) = 3
 alg_adaptive_order(alg::Rosenbrock32) = 2
 
-alg_adaptive_order(alg::RadauIIA3) = 1
-alg_adaptive_order(alg::RadauIIA5) = 3
-
 # this is actually incorrect and is purposefully decreased as this tends
 # to track the real error much better
 # this is actually incorrect and is purposefully decreased as this tends
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 1b921906bb..238926aa4c 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -554,168 +554,6 @@ struct CayleyEuler <: OrdinaryDiffEqAlgorithm end
 
 ################################################################################
 
-# FIRK Methods
-
-"""
-@article{hairer1999stiff,
-title={Stiff differential equations solved by Radau methods},
-author={Hairer, Ernst and Wanner, Gerhard},
-journal={Journal of Computational and Applied Mathematics},
-volume={111},
-number={1-2},
-pages={93--111},
-year={1999},
-publisher={Elsevier}
-}
-
-RadauIIA3: Fully-Implicit Runge-Kutta Method
-An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
-"""
-struct RadauIIA3{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
-       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
-    linsolve::F
-    precs::P
-    extrapolant::Symbol
-    κ::Tol
-    maxiters::Int
-    fast_convergence_cutoff::C1
-    new_W_γdt_cutoff::C2
-    controller::Symbol
-    step_limiter!::StepLimiter
-end
-
-function RadauIIA3(; chunk_size = Val{0}(), autodiff = Val{true}(),
-        standardtag = Val{true}(), concrete_jac = nothing,
-        diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS,
-        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
-        new_W_γdt_cutoff = 1 // 5,
-        controller = :Predictive, κ = nothing, maxiters = 10,
-        step_limiter! = trivial_limiter!)
-    RadauIIA3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(κ), typeof(fast_convergence_cutoff),
-        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
-        precs,
-        extrapolant,
-        κ,
-        maxiters,
-        fast_convergence_cutoff,
-        new_W_γdt_cutoff,
-        controller,
-        step_limiter!)
-end
-
-TruncatedStacktraces.@truncate_stacktrace RadauIIA3
-
-"""
-@article{hairer1999stiff,
-title={Stiff differential equations solved by Radau methods},
-author={Hairer, Ernst and Wanner, Gerhard},
-journal={Journal of Computational and Applied Mathematics},
-volume={111},
-number={1-2},
-pages={93--111},
-year={1999},
-publisher={Elsevier}
-}
-
-RadauIIA5: Fully-Implicit Runge-Kutta Method
-An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
-"""
-struct RadauIIA5{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
-       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
-    linsolve::F
-    precs::P
-    smooth_est::Bool
-    extrapolant::Symbol
-    κ::Tol
-    maxiters::Int
-    fast_convergence_cutoff::C1
-    new_W_γdt_cutoff::C2
-    controller::Symbol
-    step_limiter!::StepLimiter
-end
-
-function RadauIIA5(; chunk_size = Val{0}(), autodiff = Val{true}(),
-        standardtag = Val{true}(), concrete_jac = nothing,
-        diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS,
-        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
-        new_W_γdt_cutoff = 1 // 5,
-        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
-        step_limiter! = trivial_limiter!)
-    RadauIIA5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(κ), typeof(fast_convergence_cutoff),
-        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
-        precs,
-        smooth_est,
-        extrapolant,
-        κ,
-        maxiters,
-        fast_convergence_cutoff,
-        new_W_γdt_cutoff,
-        controller,
-        step_limiter!)
-end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA5
-
-"""
-@article{hairer1999stiff,
-title={Stiff differential equations solved by Radau methods},
-author={Hairer, Ernst and Wanner, Gerhard},
-journal={Journal of Computational and Applied Mathematics},
-volume={111},
-number={1-2},
-pages={93--111},
-year={1999},
-publisher={Elsevier}
-}
-
-RadauIIA7: Fully-Implicit Runge-Kutta Method
-An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
-"""
-struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
-       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
-    linsolve::F
-    precs::P
-    smooth_est::Bool
-    extrapolant::Symbol
-    κ::Tol
-    maxiters::Int
-    fast_convergence_cutoff::C1
-    new_W_γdt_cutoff::C2
-    controller::Symbol
-    step_limiter!::StepLimiter
-end
-
-function RadauIIA7(; chunk_size = Val{0}(), autodiff = Val{true}(),
-        standardtag = Val{true}(), concrete_jac = nothing,
-        diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS,
-        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
-        new_W_γdt_cutoff = 1 // 5,
-        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
-        step_limiter! = trivial_limiter!)
-    RadauIIA7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(κ), typeof(fast_convergence_cutoff),
-        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
-        precs,
-        smooth_est,
-        extrapolant,
-        κ,
-        maxiters,
-        fast_convergence_cutoff,
-        new_W_γdt_cutoff,
-        controller,
-        step_limiter!)
-end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA7
-
-################################################################################
-
 ################################################################################
 
 # Rosenbrock Methods
diff --git a/src/integrators/controllers.jl b/src/integrators/controllers.jl
index 69169566c3..10af114a40 100644
--- a/src/integrators/controllers.jl
+++ b/src/integrators/controllers.jl
@@ -393,59 +393,8 @@ else
 end
 ```
 """
-struct PredictiveController <: AbstractController
-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 alg isa Union{RadauIIA3, RadauIIA5, RadauIIA7}
-                @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 = DiffEqBase.fastpow(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) *
-               DiffEqBase.fastpow((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
+struct PredictiveController <: AbstractController
 end
 
 # Dummy controller without any method implementations.
diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl
index 58624905ac..571841f981 100644
--- a/src/integrators/integrator_interface.jl
+++ b/src/integrators/integrator_interface.jl
@@ -125,11 +125,6 @@ end
         cache::OrdinaryDiffEqMutableCache)
     (cache.tmp,)
 end
-@inline function DiffEqBase.get_tmp_cache(
-        integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA7},
-        cache::OrdinaryDiffEqMutableCache)
-    (cache.tmp, cache.atmp)
-end
 @inline function DiffEqBase.get_tmp_cache(integrator,
         alg::OrdinaryDiffEqNewtonAdaptiveAlgorithm,
         cache::OrdinaryDiffEqMutableCache)
diff --git a/test/runtests.jl b/test/runtests.jl
index 2bce643dd9..3327354672 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -77,6 +77,12 @@ function activate_sdirk()
     Pkg.instantiate()
 end
 
+function activate_firk()
+    Pkg.activate("../lib/OrdinaryDiffEqFIRK")
+    Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
+    Pkg.instantiate()
+end
+
 #Start Test Script
 
 @time begin
@@ -208,12 +214,15 @@ end
         @time @safetestset "Linear Methods Tests" include("algconvergence/linear_method_tests.jl")
         @time @safetestset "Split Methods Tests" include("algconvergence/split_methods_tests.jl")
         @time @safetestset "Rosenbrock Tests" include("algconvergence/ode_rosenbrock_tests.jl")
-        @time @safetestset "FIRK Tests" include("algconvergence/ode_firk_tests.jl")
         @time @safetestset "Linear-Nonlinear Methods Tests" include("algconvergence/linear_nonlinear_convergence_tests.jl")
         @time @safetestset "Linear-Nonlinear Krylov Methods Tests" include("algconvergence/linear_nonlinear_krylov_tests.jl")
         @time @safetestset "Quadruple precision Runge-Kutta Tests" include("algconvergence/ode_quadruple_precision_tests.jl")
     end
 
+    if !is_APPVEYOR && GROUP == "FIRK"
+        @time @safetestset "FIRK Tests" include("../lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl")
+    end
+
     if !is_APPVEYOR && GROUP == "Symplectic"
         @time @safetestset "Symplectic Tests" include("../lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl")
     end