diff --git a/Project.toml b/Project.toml
index 2a403342..ef544294 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,35 +4,19 @@ version = "5.12.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BoundaryValueDiffEqAscher = "7227322d-7511-4e07-9247-ad6ff830280e"
 BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
 BoundaryValueDiffEqFIRK = "85d9eb09-370e-4000-bb32-543851f73618"
 BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
 BoundaryValueDiffEqMIRKN = "9255f1d6-53bf-473e-b6bd-23f1ff009da4"
 BoundaryValueDiffEqShooting = "ed55bfe0-3725-4db6-871e-a1dc9f42a757"
-ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
-NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
-PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-Preferences = "21216c6a-2e73-6563-6e65-726566657250"
-RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 
 [weakdeps]
 ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c"
@@ -41,45 +25,29 @@ ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c"
 BoundaryValueDiffEqODEInterfaceExt = "ODEInterface"
 
 [compat]
-ADTypes = "1.9"
-Adapt = "4.1.1"
 Aqua = "0.8.7"
 ArrayInterface = "7.17"
-BandedMatrices = "1.7.5"
 BoundaryValueDiffEqAscher = "1"
 BoundaryValueDiffEqCore = "1.1"
 BoundaryValueDiffEqFIRK = "1.1"
 BoundaryValueDiffEqMIRK = "1.1"
 BoundaryValueDiffEqMIRKN = "1"
 BoundaryValueDiffEqShooting = "1.1"
-ConcreteStructs = "0.2.3"
 DiffEqBase = "6.158.3"
 DiffEqDevTools = "2.44"
-FastAlmostBandedMatrices = "0.1.4"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.38"
 Hwloc = "3"
 InteractiveUtils = "<0.0.1, 1"
 JET = "0.9.12"
-LineSearch = "0.1.4"
 LinearAlgebra = "1.10"
-LinearSolve = "2.36.2"
-Logging = "1.10"
-NonlinearSolveFirstOrder = "1"
 ODEInterface = "0.5"
 OrdinaryDiffEq = "6.90.1"
 Pkg = "1.10.0"
-PreallocationTools = "0.4.24"
-PrecompileTools = "1.2"
-Preferences = "1.4"
 Random = "1.10"
 ReTestItems = "1.23.1"
-RecursiveArrayTools = "3.27.0"
 Reexport = "1.2"
-SciMLBase = "2.60.0"
-Setfield = "1.1.1"
-SparseArrays = "1.10"
-SparseDiffTools = "2.23"
+SciMLBase = "2.64.0"
 StaticArrays = "1.9.8"
 Test = "1.10"
 julia = "1.10"
diff --git a/ext/BoundaryValueDiffEqODEInterfaceExt.jl b/ext/BoundaryValueDiffEqODEInterfaceExt.jl
index 617f7241..500240bb 100644
--- a/ext/BoundaryValueDiffEqODEInterfaceExt.jl
+++ b/ext/BoundaryValueDiffEqODEInterfaceExt.jl
@@ -1,24 +1,21 @@
 module BoundaryValueDiffEqODEInterfaceExt
 
-using BoundaryValueDiffEqCore, SciMLBase, ODEInterface, RecursiveArrayTools,
-      ConcreteStructs, Setfield, PreallocationTools
-
 using BoundaryValueDiffEq: BVPM2, BVPSOL, COLNEW
-import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, __extract_u0,
-                                __initial_guess_length, __extract_mesh,
-                                __flatten_initial_guess, __get_bcresid_prototype,
-                                __has_initial_guess, __initial_guess
-import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
-import ODEInterface: OptionsODE, OPT_ATOL, OPT_RTOL, OPT_METHODCHOICE, OPT_DIAGNOSTICOUTPUT,
-                     OPT_ERRORCONTROL, OPT_SINGULARTERM, OPT_MAXSTEPS, OPT_BVPCLASS,
-                     OPT_SOLMETHOD, OPT_RHS_CALLMODE, OPT_COLLOCATIONPTS, OPT_ADDGRIDPOINTS,
-                     OPT_MAXSUBINTERVALS, RHS_CALL_INSITU, evalSolution
-import ODEInterface: Bvpm2, bvpm2_init, bvpm2_solve, bvpm2_destroy, bvpm2_get_x
-import ODEInterface: bvpsol
-import ODEInterface: colnew
-
-import FastClosures: @closure
-import ForwardDiff
+using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, __extract_u0,
+                               __initial_guess_length, __extract_mesh,
+                               __flatten_initial_guess, __get_bcresid_prototype,
+                               __has_initial_guess, __initial_guess
+using SciMLBase: SciMLBase, BVProblem, TwoPointBVProblem
+using ODEInterface: OptionsODE, OPT_ATOL, OPT_RTOL, OPT_METHODCHOICE, OPT_DIAGNOSTICOUTPUT,
+                    OPT_ERRORCONTROL, OPT_SINGULARTERM, OPT_MAXSTEPS, OPT_BVPCLASS,
+                    OPT_SOLMETHOD, OPT_RHS_CALLMODE, OPT_COLLOCATIONPTS, OPT_ADDGRIDPOINTS,
+                    OPT_MAXSUBINTERVALS, RHS_CALL_INSITU, evalSolution
+using ODEInterface: Bvpm2, bvpm2_init, bvpm2_solve, bvpm2_destroy, bvpm2_get_x
+using ODEInterface: bvpsol
+using ODEInterface: colnew
+
+using FastClosures: @closure
+using ForwardDiff: ForwardDiff
 
 #------
 # BVPM2
diff --git a/lib/BoundaryValueDiffEqAscher/Project.toml b/lib/BoundaryValueDiffEqAscher/Project.toml
index 9a42e5ef..68490560 100644
--- a/lib/BoundaryValueDiffEqAscher/Project.toml
+++ b/lib/BoundaryValueDiffEqAscher/Project.toml
@@ -1,7 +1,7 @@
 name = "BoundaryValueDiffEqAscher"
 uuid = "7227322d-7511-4e07-9247-ad6ff830280e"
 authors = ["Qingyu Qu <erikqqy123@gmail.com>"]
-version = "1.1.0"
+version = "1.2.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -12,11 +12,11 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
@@ -25,7 +25,8 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 
 [compat]
 ADTypes = "1.9"
@@ -37,6 +38,7 @@ BoundaryValueDiffEqCore = "1.1.0"
 ConcreteStructs = "0.2.3"
 DiffEqBase = "6.158.3"
 DiffEqDevTools = "2.44"
+DifferentiationInterface = "0.6.24"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.38"
 Hwloc = "3"
@@ -44,7 +46,6 @@ InteractiveUtils = "<0.0.1, 1"
 JET = "0.9.12"
 LinearAlgebra = "1.10"
 LinearSolve = "2.36.2"
-Logging = "1.10"
 PreallocationTools = "0.4.24"
 PrecompileTools = "1.2"
 Preferences = "1.4"
@@ -55,7 +56,8 @@ Reexport = "1.2"
 SciMLBase = "2.59.1"
 Setfield = "1.1.1"
 SparseArrays = "1.10"
-SparseDiffTools = "2.23"
+SparseConnectivityTracer = "0.6.9"
+SparseMatrixColorings = "0.4.10"
 StaticArrays = "1.9.8"
 Test = "1.10"
 julia = "1.10"
diff --git a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
index 4284450f..31226a29 100644
--- a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
@@ -1,28 +1,31 @@
 module BoundaryValueDiffEqAscher
 
-using ADTypes
-using AlmostBlockDiagonals
-using BoundaryValueDiffEqCore
-using ConcreteStructs
-using FastClosures
-using ForwardDiff
+using ADTypes: ADTypes, AutoSparse, AutoForwardDiff
+using AlmostBlockDiagonals: AlmostBlockDiagonals, IntermediateAlmostBlockDiagonal
+using BoundaryValueDiffEqCore: BVPJacobianAlgorithm, __extract_problem_details,
+                               concrete_jacobian_algorithm, __Fix3,
+                               __concrete_nonlinearsolve_algorithm,
+                               __internal_nlsolve_problem, BoundaryValueDiffEqAlgorithm,
+                               __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
+                               __extract_mesh, get_dense_ad
+using ConcreteStructs: @concrete
+using DiffEqBase: DiffEqBase
+using DifferentiationInterface: DifferentiationInterface, Constant
+using FastClosures: @closure
+using ForwardDiff: ForwardDiff, Dual
 using LinearAlgebra
-using PreallocationTools
-using RecursiveArrayTools
-using Reexport
-using SciMLBase
-using Setfield
+using PreallocationTools: PreallocationTools, DiffCache
+using RecursiveArrayTools: VectorOfArray, recursivecopy
+using Reexport: @reexport
+using SciMLBase: SciMLBase, AbstractDiffEqInterpolation, StandardBVProblem, __solve,
+                 _unwrap_val
+using Setfield: @set!
+using SparseConnectivityTracer: SparseConnectivityTracer
+using SparseMatrixColorings: SparseMatrixColorings, GreedyColoringAlgorithm, LargestFirst
 
-import BoundaryValueDiffEqCore: BVPJacobianAlgorithm, __extract_problem_details,
-                                concrete_jacobian_algorithm, __Fix3,
-                                __concrete_nonlinearsolve_algorithm,
-                                BoundaryValueDiffEqAlgorithm, __sparse_jacobian_cache,
-                                __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
-                                __extract_mesh
+const DI = DifferentiationInterface
 
-import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
-
-@reexport using ADTypes, DiffEqBase, BoundaryValueDiffEqCore, SparseDiffTools, SciMLBase
+@reexport using ADTypes, BoundaryValueDiffEqCore, SciMLBase
 
 include("types.jl")
 include("utils.jl")
diff --git a/lib/BoundaryValueDiffEqAscher/src/ascher.jl b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
index 20b924b1..94d45315 100644
--- a/lib/BoundaryValueDiffEqAscher/src/ascher.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
@@ -136,7 +136,7 @@ function SciMLBase.__init(prob::BVProblem, alg::AbstractAscher; dt = 0.0,
     end
 
     if prob.f.bcjac === nothing
-        bcjac = construct_bc_jac(prob, bcresid_prototype, prob.problem_type)
+        bcjac = construct_bcjac(prob, bcresid_prototype)
     else
         bcjac = prob.f.bcjac
     end
@@ -315,31 +315,48 @@ function __construct_nlproblem(cache::AscherCache{iip, T}) where {iip, T}
     else
         @closure (z, p) -> @views Φ(cache, z, pt)
     end
+
     lz = reduce(vcat, cache.z)
-    sd = alg.jac_alg.diffmode isa AutoSparse ? SymbolicsSparsityDetection() :
-         NoSparsityDetection()
-    ad = alg.jac_alg.diffmode
-    lossₚ = (iip ? __Fix3 : Base.Fix2)(loss, cache.p)
-    jac_cache = __sparse_jacobian_cache(Val(iip), ad, sd, lossₚ, lz, lz)
-    jac_prototype = init_jacobian(jac_cache)
+    resid_prototype = zero(lz)
+    diffmode = if alg.jac_alg.diffmode isa AutoSparse
+        AutoSparse(get_dense_ad(alg.jac_alg.diffmode);
+            sparsity_detector = SparseConnectivityTracer.TracerSparsityDetector(),
+            coloring_algorithm = GreedyColoringAlgorithm(LargestFirst()))
+    else
+        alg.jac_alg.diffmode
+    end
+
+    jac_cache = if iip
+        DI.prepare_jacobian(loss, resid_prototype, diffmode, lz, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss, diffmode, lz, Constant(cache.p))
+    end
+
+    jac_prototype = if iip
+        DI.jacobian(loss, resid_prototype, jac_cache, diffmode, lz, Constant(cache.p))
+    else
+        DI.jacobian(loss, jac_cache, diffmode, lz, Constant(cache.p))
+    end
+
     jac = if iip
-        @closure (J, u, p) -> __ascher_mpoint_jacobian!(J, u, ad, jac_cache, lossₚ, lz)
+        @closure (J, u, p) -> __ascher_mpoint_jacobian!(
+            J, u, diffmode, jac_cache, loss, lz, cache.p)
     else
-        @closure (u, p) -> __ascher_mpoint_jacobian(jac_prototype, u, ad, jac_cache, lossₚ)
+        @closure (u, p) -> __ascher_mpoint_jacobian(
+            jac_prototype, u, diffmode, jac_cache, loss, cache.p)
     end
-    resid_prototype = zero(lz)
-    _nlf = NonlinearFunction{iip}(
+
+    nlf = NonlinearFunction{iip}(
         loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    nlprob::NonlinearProblem = NonlinearProblem(_nlf, lz, cache.p)
-    return nlprob
+    return __internal_nlsolve_problem(cache.prob, similar(lz), lz, nlf, lz, cache.p)
 end
 
-function __ascher_mpoint_jacobian!(J, x, diffmode, diffcache, loss, resid)
-    sparse_jacobian!(J, diffmode, diffcache, loss, resid, x)
+function __ascher_mpoint_jacobian!(J, x, diffmode, diffcache, loss, resid, p)
+    DI.jacobian!(loss, resid, J, diffcache, diffmode, x, Constant(p))
     return nothing
 end
-function __ascher_mpoint_jacobian(J, x, diffmode, diffcache, loss)
-    sparse_jacobian!(J, diffmode, diffcache, loss, x)
+function __ascher_mpoint_jacobian(J, x, diffmode, diffcache, loss, p)
+    DI.jacobian!(loss, J, diffcache, diffmode, x, Constant(p))
     return J
 end
 
diff --git a/lib/BoundaryValueDiffEqAscher/src/collocation.jl b/lib/BoundaryValueDiffEqAscher/src/collocation.jl
index cd9c0e42..592606dd 100644
--- a/lib/BoundaryValueDiffEqAscher/src/collocation.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/collocation.jl
@@ -313,9 +313,10 @@ function Φ!(cache::AscherCache{iip, T}, z, res, pt::TwoPointBVProblem) where {i
     copyto!(cache.dmz, dmz)
 end
 
-@inline __get_value(z::Vector{<:AbstractArray}) = eltype(first(z)) <: ForwardDiff.Dual ?
-                                                  [map(x -> x.value, a) for a in z] : z
-@inline __get_value(z) = isa(z, ForwardDiff.Dual) ? z.value : z
+@inline __get_value(x) = x
+@inline __get_value(x::Dual) = ForwardDiff.value(x)
+@inline __get_value(x::AbstractArray{<:Dual}) = map(ForwardDiff.value, x)
+@inline __get_value(x::AbstractArray{<:AbstractArray{<:Dual}}) = map(__get_value, x)
 
 function Φ(cache::AscherCache{iip, T}, z, pt::StandardBVProblem) where {iip, T}
     (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, residual, zval, yval, gval, delz, dmz, deldmz, g, w, v, dmzo, ipvtg, ipvtw, TU) = cache
@@ -374,7 +375,7 @@ function Φ(cache::AscherCache{iip, T}, z, pt::StandardBVProblem) where {iip, T}
 
         @views gblock!(cache, h, g[i], izeta, w[i], v[i])
 
-        if i >= n
+        if i == n
             izsave = izeta
             # build equation for a side condition.
             # other nonlinear case
@@ -739,7 +740,6 @@ end
 function vwblok(cache::AscherCache, xcol, hrho, jj, wi, vi, ipvtw, zyval, df, acol, dmzo)
     (; jac, k, p, ncomp, ny) = cache
     ncy = ncomp + ny
-    kdy = k * ncy
     # initialize wi
     i0 = (jj - 1) * ncy
     for id in (i0 + 1):(i0 + ncomp)
diff --git a/lib/BoundaryValueDiffEqAscher/src/utils.jl b/lib/BoundaryValueDiffEqAscher/src/utils.jl
index 01142a8b..bf9e6b59 100644
--- a/lib/BoundaryValueDiffEqAscher/src/utils.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/utils.jl
@@ -174,7 +174,9 @@ end
     return nothing
 end
 
-@inline function construct_bc_jac(prob::BVProblem, _, pt::StandardBVProblem)
+@inline construct_bcjac(prob, bcresid_prototype) = construct_bcjac(
+    prob, bcresid_prototype, prob.problem_type)
+@inline function construct_bcjac(prob::BVProblem, _, pt::StandardBVProblem)
     if isinplace(prob)
         bcjac = (df, u, p, t) -> begin
             _du = similar(u)
@@ -184,17 +186,15 @@ end
             return
         end
     else
-        bcjac = (df, u, p, t) -> begin
-            _du = prob.f.bc(u, p, t)
+        bcjac = (u, p, t) -> begin
             _f = @closure (du, u) -> (du .= prob.f.bc(u, p, t))
-            ForwardDiff.jacobian!(df, _f, _du, u)
-            return
+            return ForwardDiff.jacobian(_f, u)
         end
     end
     return bcjac
 end
 
-@inline function construct_bc_jac(prob::BVProblem, bcresid_prototype, pt::TwoPointBVProblem)
+@inline function construct_bcjac(prob::BVProblem, bcresid_prototype, pt::TwoPointBVProblem)
     if isinplace(prob)
         bcjac = (df, u, p) -> begin
             _du = similar(u)
diff --git a/lib/BoundaryValueDiffEqCore/Project.toml b/lib/BoundaryValueDiffEqCore/Project.toml
index 414bc306..cb079ae3 100644
--- a/lib/BoundaryValueDiffEqCore/Project.toml
+++ b/lib/BoundaryValueDiffEqCore/Project.toml
@@ -21,7 +21,6 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 
 [compat]
 ADTypes = "1.9"
@@ -43,7 +42,6 @@ Reexport = "1.2"
 SciMLBase = "2.59.1"
 Setfield = "1"
 SparseArrays = "1.10"
-SparseDiffTools = "2.23"
 Test = "1.10"
 julia = "1.10"
 
diff --git a/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl b/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl
index bf68de0a..ecedd838 100644
--- a/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl
+++ b/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl
@@ -1,36 +1,35 @@
 module BoundaryValueDiffEqCore
 
-using ADTypes, Adapt, ArrayInterface, ForwardDiff, LinearAlgebra, LineSearch,
-      NonlinearSolveFirstOrder, RecursiveArrayTools, Reexport, SciMLBase, Setfield,
-      SparseDiffTools
-
-using PreallocationTools: PreallocationTools, DiffCache
-
-# Special Matrix Types
-using SparseArrays
-
-import ADTypes: AbstractADType
-import ArrayInterface: matrix_colors, parameterless_type, fast_scalar_indexing
-import ConcreteStructs: @concrete
-import DiffEqBase: solve
-import ForwardDiff: ForwardDiff, pickchunksize, Dual
-import Logging
+using Adapt: adapt
+using ADTypes: ADTypes, AbstractADType, AutoSparse, AutoForwardDiff, AutoFiniteDiff,
+               NoSparsityDetector, KnownJacobianSparsityDetector
+using ArrayInterface: matrix_colors, parameterless_type, fast_scalar_indexing
+using ConcreteStructs: @concrete
+using DiffEqBase: DiffEqBase, solve
+using ForwardDiff: ForwardDiff, pickchunksize
+using Logging
 using NonlinearSolveFirstOrder: NonlinearSolvePolyAlgorithm
-import LineSearch: BackTracking
-import RecursiveArrayTools: VectorOfArray, DiffEqArray
-import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
+using LinearAlgebra
+using LineSearch: BackTracking
+using PreallocationTools: PreallocationTools, DiffCache
+using RecursiveArrayTools: AbstractVectorOfArray, VectorOfArray, DiffEqArray
+using Reexport: @reexport
+using SciMLBase: SciMLBase, AbstractBVProblem, AbstractDiffEqInterpolation,
+                 StandardBVProblem, StandardSecondOrderBVProblem, __solve, _unwrap_val
+using Setfield: @set!, @set
+using SparseArrays: sparse
 
-@reexport using ADTypes, NonlinearSolveFirstOrder, SparseDiffTools, SciMLBase
+@reexport using NonlinearSolveFirstOrder, SciMLBase
 
 include("types.jl")
 include("utils.jl")
 include("algorithms.jl")
 include("alg_utils.jl")
 include("default_nlsolve.jl")
-include("sparse_jacobians.jl")
 include("misc_utils.jl")
 
-function __solve(prob::BVProblem, alg::BoundaryValueDiffEqAlgorithm, args...; kwargs...)
+function SciMLBase.__solve(
+        prob::AbstractBVProblem, alg::BoundaryValueDiffEqAlgorithm, args...; kwargs...)
     cache = init(prob, alg, args...; kwargs...)
     return solve!(cache)
 end
diff --git a/lib/BoundaryValueDiffEqCore/src/sparse_jacobians.jl b/lib/BoundaryValueDiffEqCore/src/sparse_jacobians.jl
deleted file mode 100644
index ac4e571e..00000000
--- a/lib/BoundaryValueDiffEqCore/src/sparse_jacobians.jl
+++ /dev/null
@@ -1,38 +0,0 @@
-# This file defines several common patterns of sparse Jacobians we see in the BVP solvers.
-function _sparse_like(I, J, x::AbstractArray, m = maximum(I), n = maximum(J))
-    I′ = adapt(parameterless_type(x), I)
-    J′ = adapt(parameterless_type(x), J)
-    V = __ones_like(x, length(I))
-    return sparse(I′, J′, V, m, n)
-end
-
-# NOTE: We don't retain the Banded Structure in non-TwoPoint BVP cases since vcat/hcat makes
-# it into a dense array. Instead we can atleast exploit sparsity!
-
-# FIXME: Fix the cases where fast_scalar_indexing is not possible
-
-# Helpers for IIP/OOP functions
-function __sparse_jacobian_cache(::Val{iip}, ad, sd, fn, fx, y) where {iip}
-    if iip
-        sparse_jacobian_cache(ad, sd, fn, fx, y)
-    else
-        sparse_jacobian_cache(ad, sd, fn, y; fx)
-    end
-end
-
-@concrete struct ColoredMatrix
-    M
-    row_colorvec
-    col_colorvec
-end
-
-Base.size(M::ColoredMatrix, args...) = size(M.M, args...)
-Base.eltype(M::ColoredMatrix) = eltype(M.M)
-
-ColoredMatrix() = ColoredMatrix(nothing, nothing, nothing)
-
-function __sparsity_detection_alg(M::ColoredMatrix)
-    return PrecomputedJacobianColorvec(;
-        jac_prototype = M.M, M.row_colorvec, M.col_colorvec)
-end
-__sparsity_detection_alg(::ColoredMatrix{Nothing}) = NoSparsityDetection()
diff --git a/lib/BoundaryValueDiffEqCore/src/types.jl b/lib/BoundaryValueDiffEqCore/src/types.jl
index b1f34007..5c4910b8 100644
--- a/lib/BoundaryValueDiffEqCore/src/types.jl
+++ b/lib/BoundaryValueDiffEqCore/src/types.jl
@@ -43,7 +43,7 @@ function BVPJacobianAlgorithm(
     if diffmode !== missing
         bc_diffmode = bc_diffmode === missing ? diffmode : bc_diffmode
         nonbc_diffmode = nonbc_diffmode === missing ? diffmode : nonbc_diffmode
-        return BVPJacobianAlgorithm(diffmode, diffmode, diffmode)
+        return BVPJacobianAlgorithm(bc_diffmode, nonbc_diffmode, diffmode)
     else
         diffmode = nothing
         bc_diffmode = bc_diffmode === missing ? nothing : bc_diffmode
diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 33afd602..f786547b 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -104,6 +104,47 @@ end
     return resid
 end
 
+function eval_bc_residual(::StandardSecondOrderBVProblem, bc::BC, y, dy, p, mesh) where {BC}
+    L = length(mesh)
+    res_bc = bc(dy, y, p, mesh)
+    return res_bc
+end
+function eval_bc_residual(
+        ::TwoPointSecondOrderBVProblem, (bca, bcb)::BC, sol, p, mesh) where {BC}
+    M = length(sol[1])
+    L = length(mesh)
+    ua = sol isa VectorOfArray ? sol[:, 1] : sol(first(t))[1:M]
+    ub = sol isa VectorOfArray ? sol[:, L] : sol(last(t))[1:M]
+    dua = sol isa VectorOfArray ? sol[:, L + 1] : sol(first(t))[(M + 1):end]
+    dub = sol isa VectorOfArray ? sol[:, end] : sol(last(t))[(M + 1):end]
+    return vcat(bca(dua, ua, p), bcb(dub, ub, p))
+end
+
+function eval_bc_residual!(resid, ::StandardBVProblem, bc!::BC, sol, p, t) where {BC}
+    bc!(resid, sol, p, t)
+end
+
+function eval_bc_residual!(
+        resid, ::StandardSecondOrderBVProblem, bc!::BC, sol, dsol, p, mesh) where {BC}
+    M = length(sol[1])
+    res_bc = vcat(resid[1], resid[2])
+    bc!(res_bc, dsol, sol, p, mesh)
+    copyto!(resid[1], res_bc[1:M])
+    copyto!(resid[2], res_bc[(M + 1):end])
+end
+
+function eval_bc_residual!(
+        resid, ::TwoPointSecondOrderBVProblem, (bca!, bcb!)::BC, sol, p, mesh) where {BC}
+    M = length(sol[1])
+    L = length(mesh)
+    ua = sol isa VectorOfArray ? sol[:, 1] : sol(first(mesh))[1:M]
+    ub = sol isa VectorOfArray ? sol[:, L] : sol(last(mesh))[1:M]
+    dua = sol isa VectorOfArray ? sol[:, L + 1] : sol(first(mesh))[(M + 1):end]
+    dub = sol isa VectorOfArray ? sol[:, end] : sol(last(mesh))[(M + 1):end]
+    bca!(resid[1], dua, ua, p)
+    bcb!(resid[2], dub, ub, p)
+end
+
 __append_similar!(::Nothing, n, _) = nothing
 
 # NOTE: We use `last` since the `first` might not conform to the same structure. For eg,
@@ -143,8 +184,7 @@ function __extract_problem_details(prob, u0::AbstractVector{<:AbstractArray}; kw
     _u0 = first(u0)
     return Val(true), eltype(_u0), length(_u0), (length(u0) - 1), _u0
 end
-function __extract_problem_details(
-        prob, u0::RecursiveArrayTools.AbstractVectorOfArray; kwargs...)
+function __extract_problem_details(prob, u0::AbstractVectorOfArray; kwargs...)
     # Problem has Initial Guess
     _u0 = first(u0.u)
     return Val(true), eltype(_u0), length(_u0), (length(u0.u) - 1), _u0
@@ -195,6 +235,20 @@ function __get_bcresid_prototype(::StandardBVProblem, prob::BVProblem, u)
     return prototype, size(prototype)
 end
 
+function __get_bcresid_prototype(::TwoPointSecondOrderBVProblem, prob::BVProblem, u)
+    prototype = if prob.f.bcresid_prototype !== nothing
+        prob.f.bcresid_prototype.x
+    else
+        first(prob.f.bc)(u, prob.p), last(prob.f.bc)(u, prob.p)
+    end
+    return prototype, size.(prototype)
+end
+function __get_bcresid_prototype(::StandardSecondOrderBVProblem, prob::BVProblem, u)
+    prototype = prob.f.bcresid_prototype !== nothing ? prob.f.bcresid_prototype :
+                __zeros_like(u)
+    return prototype, size(prototype)
+end
+
 @inline function __similar(x, args...)
     y = similar(x, args...)
     return zero(y)
@@ -235,6 +289,32 @@ end
 __vec_bc(sol, p, t, bc, u_size) = vec(bc(sol, p, t))
 __vec_bc(sol, p, bc, u_size) = vec(bc(reshape(sol, u_size), p))
 
+# Restructure Non-Vector Inputs
+function __vec_f!(ddu, du, u, p, t, f!, u_size)
+    f!(reshape(ddu, u_size), reshape(du, u_size), reshape(u, u_size), p, t)
+    return nothing
+end
+
+__vec_f(du, u, p, t, f, u_size) = vec(f(reshape(du, u_size), reshape(u, u_size), p, t))
+
+function __vec_so_bc!(resid, dsol, sol, p, t, bc!, resid_size, u_size)
+    bc!(reshape(resid, resid_size), __restructure_sol(dsol, u_size),
+        __restructure_sol(sol, u_size), p, t)
+    return nothing
+end
+
+function __vec_so_bc!(resid, dsol, sol, p, bc!, resid_size, u_size)
+    bc!(reshape(resid, resid_size), reshape(dsol, u_size), reshape(sol, u_size), p)
+    return nothing
+end
+
+function __vec_so_bc(dsol, sol, p, t, bc, u_size)
+    vec(bc(__restructure_sol(dsol, u_size), __restructure_sol(sol, u_size), p, t))
+end
+function __vec_so_bc(dsol, sol, p, bc, u_size)
+    vec(bc(reshape(dsol, u_size), reshape(sol, u_size), p))
+end
+
 @inline __get_non_sparse_ad(ad::AbstractADType) = ad
 @inline __get_non_sparse_ad(ad::AutoSparse) = ADTypes.dense_ad(ad)
 
@@ -276,6 +356,12 @@ end
     end
 end
 
+@inline function __internal_nlsolve_problem(
+        ::SecondOrderBVProblem{uType, tType, iip, nlls}, resid_prototype,
+        u0, args...; kwargs...) where {uType, tType, iip, nlls}
+    return NonlinearProblem(args...; kwargs...)
+end
+
 # Handling Initial Guesses
 """
     __extract_u0(u₀, t₀)
@@ -356,6 +442,10 @@ end
 @inline function __initial_guess_on_mesh(u₀::F, mesh, p) where {F}
     return VectorOfArray([vec(__initial_guess(u₀, p, t)) for t in mesh])
 end
+@inline function __initial_guess_on_mesh(
+        prob::SecondOrderBVProblem, u₀::AbstractArray, Nig, p, alias_u0::Bool)
+    return VectorOfArray([copy(vec(u₀)) for _ in 1:(2 * (Nig + 1))])
+end
 
 # Construct BVP Solution
 function __build_solution(prob::BVProblem, odesol, nlsol)
@@ -370,3 +460,18 @@ end
 end
 
 @inline (f::__Fix3{F})(a, b) where {F} = f.f(a, b, f.x)
+
+# convert every vector of vector to AbstractVectorOfArray, especially if them come from get_tmp of PreallocationTools.jl
+
+get_dense_ad(::Nothing) = nothing
+get_dense_ad(ad) = ad
+get_dense_ad(ad::AutoSparse) = ADTypes.dense_ad(ad)
+
+# traits for forward or reverse mode AutoForwardDiff
+
+function _sparse_like(I, J, x::AbstractArray, m = maximum(I), n = maximum(J))
+    I′ = adapt(parameterless_type(x), I)
+    J′ = adapt(parameterless_type(x), J)
+    V = __ones_like(x, length(I))
+    return sparse(I′, J′, V, m, n)
+end
diff --git a/lib/BoundaryValueDiffEqFIRK/Project.toml b/lib/BoundaryValueDiffEqFIRK/Project.toml
index ce9571fc..8262e0eb 100644
--- a/lib/BoundaryValueDiffEqFIRK/Project.toml
+++ b/lib/BoundaryValueDiffEqFIRK/Project.toml
@@ -1,21 +1,20 @@
 name = "BoundaryValueDiffEqFIRK"
 uuid = "85d9eb09-370e-4000-bb32-543851f73618"
-version = "1.2.0"
+version = "1.3.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
@@ -24,11 +23,11 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 
 [compat]
 ADTypes = "1.9"
-Adapt = "4.1.1"
 Aqua = "0.8.7"
 ArrayInterface = "7.16"
 BandedMatrices = "1.7.5"
@@ -36,6 +35,7 @@ BoundaryValueDiffEqCore = "1.1.0"
 ConcreteStructs = "0.2.3"
 DiffEqBase = "6.158.3"
 DiffEqDevTools = "2.44"
+DifferentiationInterface = "0.6.24"
 FastAlmostBandedMatrices = "0.1.4"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.38"
@@ -44,7 +44,6 @@ InteractiveUtils = "<0.0.1, 1"
 JET = "0.9.12"
 LinearAlgebra = "1.10"
 LinearSolve = "2.36.2"
-Logging = "1.10"
 OrdinaryDiffEq = "6.90.1"
 PreallocationTools = "0.4.24"
 PrecompileTools = "1.2"
@@ -53,10 +52,11 @@ Random = "1.10"
 ReTestItems = "1.23.1"
 RecursiveArrayTools = "3.27.0"
 Reexport = "1.2"
-SciMLBase = "2.59.1"
+SciMLBase = "2.64.0"
 Setfield = "1.1.1"
 SparseArrays = "1.10"
-SparseDiffTools = "2.23"
+SparseConnectivityTracer = "0.6.9"
+SparseMatrixColorings = "0.4.10"
 StaticArrays = "1.9.8"
 Test = "1.10"
 julia = "1.10"
diff --git a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
index 860a3af9..11e8b35c 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
@@ -1,45 +1,51 @@
 module BoundaryValueDiffEqFIRK
 
-import PrecompileTools: @compile_workload, @setup_workload
-
-using ADTypes, Adapt, ArrayInterface, BoundaryValueDiffEqCore, DiffEqBase, ForwardDiff,
-      LinearAlgebra, Preferences, RecursiveArrayTools, Reexport, SciMLBase, Setfield,
-      SparseDiffTools
-
+using ADTypes: ADTypes, AutoSparse, AutoForwardDiff
+using ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
+using BandedMatrices: BandedMatrix, Ones
+using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
+                               recursive_flatten, recursive_flatten!, recursive_unflatten!,
+                               __concrete_nonlinearsolve_algorithm, diff!,
+                               __FastShortcutBVPCompatibleNonlinearPolyalg,
+                               __FastShortcutBVPCompatibleNLLSPolyalg, EvalSol,
+                               concrete_jacobian_algorithm, eval_bc_residual,
+                               eval_bc_residual!, get_tmp, __maybe_matmul!,
+                               __append_similar!, __extract_problem_details,
+                               __initial_guess, __maybe_allocate_diffcache,
+                               __restructure_sol, __get_bcresid_prototype, __similar, __vec,
+                               __vec_f, __vec_f!, __vec_bc, __vec_bc!,
+                               recursive_flatten_twopoint!, __internal_nlsolve_problem,
+                               MaybeDiffCache, __extract_mesh, __extract_u0,
+                               __has_initial_guess, __initial_guess_length,
+                               __initial_guess_on_mesh, __flatten_initial_guess,
+                               __build_solution, __Fix3, get_dense_ad, _sparse_like
+
+using ConcreteStructs: @concrete
+using DiffEqBase: DiffEqBase
+using DifferentiationInterface: DifferentiationInterface, Constant
+using FastAlmostBandedMatrices: AlmostBandedMatrix, fillpart, exclusive_bandpart,
+                                finish_part_setindex!
+using FastClosures: @closure
+using ForwardDiff: ForwardDiff, pickchunksize, Dual
+using LinearAlgebra
+using RecursiveArrayTools: AbstractVectorOfArray, AbstractVectorOfArray, DiffEqArray,
+                           VectorOfArray, recursivecopy, recursivefill!
+using Reexport: @reexport
 using PreallocationTools: PreallocationTools, DiffCache
-
-# Special Matrix Types
-using BandedMatrices, FastAlmostBandedMatrices, SparseArrays
-
-import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
-                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                                __concrete_nonlinearsolve_algorithm, diff!,
-                                __FastShortcutBVPCompatibleNonlinearPolyalg,
-                                __FastShortcutBVPCompatibleNLLSPolyalg, EvalSol,
-                                concrete_jacobian_algorithm, eval_bc_residual,
-                                eval_bc_residual!, get_tmp, __maybe_matmul!,
-                                __append_similar!, __extract_problem_details,
-                                __initial_guess, __maybe_allocate_diffcache,
-                                __restructure_sol, __get_bcresid_prototype, __similar,
-                                __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
-                                recursive_flatten_twopoint!, __internal_nlsolve_problem,
-                                MaybeDiffCache, __extract_mesh, __extract_u0,
-                                __has_initial_guess, __initial_guess_length,
-                                __initial_guess_on_mesh, __flatten_initial_guess,
-                                __build_solution, __Fix3, __sparse_jacobian_cache,
-                                __sparsity_detection_alg, _sparse_like, ColoredMatrix
-
-import ADTypes: AbstractADType
-import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
-import ConcreteStructs: @concrete
-import DiffEqBase: solve
-import FastClosures: @closure
-import ForwardDiff: ForwardDiff, pickchunksize, Dual
-import Logging
-import RecursiveArrayTools: ArrayPartition, DiffEqArray
-import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
-
-@reexport using ADTypes, DiffEqBase, BoundaryValueDiffEqCore, SparseDiffTools, SciMLBase
+using PrecompileTools: @compile_workload, @setup_workload
+using Preferences: Preferences
+using SciMLBase: SciMLBase, AbstractDiffEqInterpolation, StandardBVProblem, __solve,
+                 _unwrap_val
+using Setfield: @set!, @set
+using SparseArrays: sparse
+using SparseConnectivityTracer: SparseConnectivityTracer
+using SparseMatrixColorings: ColoringProblem, GreedyColoringAlgorithm, sparsity_pattern,
+                             ConstantColoringAlgorithm, row_colors, column_colors, coloring,
+                             LargestFirst
+
+const DI = DifferentiationInterface
+
+@reexport using ADTypes, BoundaryValueDiffEqCore, SciMLBase
 
 include("types.jl")
 include("utils.jl")
diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 404a9dae..52fd073a 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -444,43 +444,66 @@ function __construct_nlproblem(
     L = length(resid_bc)
     resid_collocation = __similar(y, cache.M * (N - 1) * (stage + 1))
 
-    loss_bcₚ = (iip ? __Fix3 : Base.Fix2)(loss_bc, cache.p)
-    loss_collocationₚ = (iip ? __Fix3 : Base.Fix2)(loss_collocation, cache.p)
+    bc_diffmode = if jac_alg.bc_diffmode isa AutoSparse
+        AutoSparse(jac_alg.bc_diffmode;
+            sparsity_detector = SparseConnectivityTracer.TracerSparsityDetector(),
+            coloring_algorithm = GreedyColoringAlgorithm())
+    else
+        jac_alg.bc_diffmode
+    end
 
-    sd_bc = jac_alg.bc_diffmode isa AutoSparse ? SymbolicsSparsityDetection() :
-            NoSparsityDetection()
-    cache_bc = __sparse_jacobian_cache(
-        Val(iip), jac_alg.bc_diffmode, sd_bc, loss_bcₚ, resid_bc, y)
+    cache_bc = if iip
+        DI.prepare_jacobian(loss_bc, resid_bc, bc_diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss_bc, bc_diffmode, y, Constant(cache.p))
+    end
 
-    sd_collocation = if jac_alg.nonbc_diffmode isa AutoSparse
+    nonbc_diffmode = if jac_alg.nonbc_diffmode isa AutoSparse
         if L < cache.M
             # For underdetermined problems we use sparse since we don't have banded qr
-            colored_matrix = __generate_sparse_jacobian_prototype(
-                cache, cache.problem_type, y, y, cache.M, N)
             J_full_band = nothing
-            __sparsity_detection_alg(ColoredMatrix(
-                sparse(colored_matrix.M), colored_matrix.row_colorvec,
-                colored_matrix.col_colorvec))
+            colored_result = __generate_sparse_jacobian_prototype(
+                cache, cache.problem_type, y, y, cache.M, N, jac_alg.nonbc_diffmode)
         else
             block_size = cache.M * (stage + 2)
             J_full_band = BandedMatrix(
                 Ones{eltype(y)}(L + cache.M * (stage + 1) * (N - 1),
                     cache.M * (stage + 1) * (N - 1) + cache.M),
                 (block_size, block_size))
-            __sparsity_detection_alg(__generate_sparse_jacobian_prototype(
-                cache, cache.problem_type, y, y, cache.M, N))
+            colored_result = __generate_sparse_jacobian_prototype(
+                cache, cache.problem_type, y, y, cache.M, N, jac_alg.nonbc_diffmode)
         end
+        is_diffmode_reverse = ADTypes.mode(jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode
+        partition = ifelse(is_diffmode_reverse, :row, :column)
+        constant_matrix_coloring = ifelse(is_diffmode_reverse, row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(jac_alg.nonbc_diffmode);
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
     else
         J_full_band = nothing
-        NoSparsityDetection()
+        jac_alg.nonbc_diffmode
     end
 
-    cache_collocation = __sparse_jacobian_cache(
-        Val(iip), jac_alg.nonbc_diffmode, sd_collocation,
-        loss_collocationₚ, resid_collocation, y)
+    cache_collocation = if iip
+        DI.prepare_jacobian(
+            loss_collocation, resid_collocation, nonbc_diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss_collocation, nonbc_diffmode, y, Constant(cache.p))
+    end
 
-    J_bc = zero(init_jacobian(cache_bc))
-    J_c = zero(init_jacobian(cache_collocation))
+    J_bc = if iip
+        DI.jacobian(loss_bc, resid_bc, cache_bc, bc_diffmode, y, Constant(cache.p))
+    else
+        DI.jacobian(loss_bc, cache_bc, bc_diffmode, y, Constant(cache.p))
+    end
+    J_c = if iip
+        DI.jacobian(loss_collocation, resid_collocation, cache_collocation,
+            nonbc_diffmode, y, Constant(cache.p))
+    else
+        DI.jacobian(
+            loss_collocation, cache_collocation, nonbc_diffmode, y, Constant(cache.p))
+    end
 
     if J_full_band === nothing
         jac_prototype = vcat(J_bc, J_c)
@@ -490,16 +513,14 @@ function __construct_nlproblem(
 
     jac = if iip
         @closure (J, u, p) -> __firk_mpoint_jacobian!(
-            J, J_c, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode, cache_bc,
-            cache_collocation, loss_bcₚ, loss_collocationₚ, resid_bc, resid_collocation, L)
+            J, J_c, u, bc_diffmode, nonbc_diffmode, cache_bc, cache_collocation,
+            loss_bc, loss_collocation, resid_bc, resid_collocation, L, cache.p)
     else
         @closure (u, p) -> __firk_mpoint_jacobian(
-            jac_prototype, J_c, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode,
-            cache_bc, cache_collocation, loss_bcₚ, loss_collocationₚ, L)
+            jac_prototype, J_c, u, bc_diffmode, nonbc_diffmode, cache_bc,
+            cache_collocation, loss_bc, loss_collocation, L, cache.p)
     end
 
-    resid_prototype = vcat(resid_bc, resid_collocation)
-
     resid_prototype = vcat(resid_bc, resid_collocation)
     nlf = NonlinearFunction{iip}(
         loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
@@ -514,8 +535,6 @@ function __construct_nlproblem(
     (; stage) = cache
     N = length(cache.mesh)
 
-    lossₚ = iip ? ((du, u) -> loss(du, u, cache.p)) : (u -> loss(u, cache.p))
-
     resid_collocation = __similar(y, cache.M * (N - 1) * (stage + 1))
 
     resid = vcat(
@@ -523,31 +542,47 @@ function __construct_nlproblem(
         @view(cache.bcresid_prototype[(prod(cache.resid_size[1]) + 1):end]))
     L = length(cache.bcresid_prototype)
 
-    sd = if jac_alg.nonbc_diffmode isa AutoSparse
+    diffmode = if jac_alg.diffmode isa AutoSparse
         block_size = cache.M * (stage + 2)
         J_full_band = BandedMatrix(
             Ones{eltype(y)}(L + cache.M * (stage + 1) * (N - 1),
                 cache.M * (stage + 1) * (N - 1) + cache.M),
             (block_size, block_size))
-        __sparsity_detection_alg(__generate_sparse_jacobian_prototype(
-            cache, cache.problem_type,
+        colored_result = __generate_sparse_jacobian_prototype(cache, cache.problem_type,
             @view(cache.bcresid_prototype[1:prod(cache.resid_size[1])]),
             @view(cache.bcresid_prototype[(prod(cache.resid_size[1]) + 1):end]),
-            cache.M, N))
+            cache.M, N, jac_alg.diffmode)
+        partition = ifelse(
+            ADTypes.mode(jac_alg.bc_diffmode) isa ADTypes.ReverseMode, :row, :column)
+        constant_matrix_coloring = ifelse(
+            ADTypes.mode(jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode,
+            row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(jac_alg.diffmode);
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
     else
-        J_full_band = nothing
-        NoSparsityDetection()
+        jac_alg.diffmode
     end
 
-    diffcache = __sparse_jacobian_cache(Val(iip), jac_alg.diffmode, sd, lossₚ, resid, y)
-    jac_prototype = zero(init_jacobian(diffcache))
+    diffcache = if iip
+        DI.prepare_jacobian(loss, resid, diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss, diffmode, y, Constant(cache.p))
+    end
+
+    jac_prototype = if iip
+        DI.jacobian(loss, resid, diffcache, diffmode, y, Constant(cache.p)) #zero(init_jacobian(diffcache))
+    else
+        DI.jacobian(loss, diffcache, diffmode, y, Constant(cache.p))
+    end
 
     jac = if iip
         @closure (J, u, p) -> __firk_2point_jacobian!(
-            J, u, jac_alg.diffmode, diffcache, lossₚ, resid)
+            J, u, jac_alg.diffmode, diffcache, loss, resid, cache.p)
     else
         @closure (u, p) -> __firk_2point_jacobian(
-            u, jac_prototype, jac_alg.diffmode, diffcache, lossₚ)
+            u, jac_prototype, jac_alg.diffmode, diffcache, loss, cache.p)
     end
 
     resid_prototype = copy(resid)
@@ -566,39 +601,66 @@ function __construct_nlproblem(
     L = length(resid_bc)
     resid_collocation = __similar(y, cache.M * (N - 1))
 
-    loss_bcₚ = (iip ? __Fix3 : Base.Fix2)(loss_bc, cache.p)
-    loss_collocationₚ = (iip ? __Fix3 : Base.Fix2)(loss_collocation, cache.p)
+    bc_diffmode = if jac_alg.bc_diffmode isa AutoSparse
+        AutoSparse(get_dense_ad(jac_alg.bc_diffmode);
+            sparsity_detector = SparseConnectivityTracer.TracerSparsityDetector(),
+            coloring_algorithm = GreedyColoringAlgorithm())
+    else
+        jac_alg.bc_diffmode
+    end
 
-    sd_bc = jac_alg.bc_diffmode isa AutoSparse ? SymbolicsSparsityDetection() :
-            NoSparsityDetection()
-    cache_bc = __sparse_jacobian_cache(
-        Val(iip), jac_alg.bc_diffmode, sd_bc, loss_bcₚ, resid_bc, y)
+    cache_bc = if iip
+        DI.prepare_jacobian(loss_bc, resid_bc, bc_diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss_bc, bc_diffmode, y, Constant(cache.p))
+    end
 
-    sd_collocation = if jac_alg.nonbc_diffmode isa AutoSparse
+    nonbc_diffmode = if jac_alg.nonbc_diffmode isa AutoSparse
         if L < cache.M
             # For underdetermined problems we use sparse since we don't have banded qr
-            colored_matrix = __generate_sparse_jacobian_prototype(
-                cache, cache.problem_type, y, y, cache.M, N)
             J_full_band = nothing
-            __sparsity_detection_alg(ColoredMatrix(
-                sparse(colored_matrix.M), colored_matrix.row_colorvec,
-                colored_matrix.col_colorvec))
+            colored_result = __generate_sparse_jacobian_prototype(
+                cache, cache.problem_type, y, y, cache.M, N, jac_alg.nonbc_diffmode)
         else
             J_full_band = BandedMatrix(Ones{eltype(y)}(L + cache.M * (N - 1), cache.M * N),
                 (L + 1, cache.M + max(cache.M - L, 0)))
-            __sparsity_detection_alg(__generate_sparse_jacobian_prototype(
-                cache, cache.problem_type, y, y, cache.M, N))
+            colored_result = __generate_sparse_jacobian_prototype(
+                cache, cache.problem_type, y, y, cache.M, N, jac_alg.nonbc_diffmode)
         end
+        partition = ifelse(
+            ADTypes.mode(jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode, :row, :column)
+        constant_matrix_coloring = ifelse(
+            ADTypes.mode(jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode,
+            row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(jac_alg.nonbc_diffmode);
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
     else
         J_full_band = nothing
-        NoSparsityDetection()
+        jac_alg.nonbc_diffmode
+    end
+
+    cache_collocation = if iip
+        DI.prepare_jacobian(
+            loss_collocation, resid_collocation, nonbc_diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss_collocation, nonbc_diffmode, y, Constant(cache.p))
+    end
+
+    J_bc = if iip
+        DI.jacobian(loss_bc, resid_bc, cache_bc, bc_diffmode, y, Constant(cache.p))
+    else
+        DI.jacobian(loss_bc, cache_bc, bc_diffmode, y, Constant(cache.p))
+    end
+    J_c = if iip
+        DI.jacobian(loss_collocation, resid_collocation, cache_collocation,
+            nonbc_diffmode, y, Constant(cache.p))
+    else
+        DI.jacobian(
+            loss_collocation, cache_collocation, nonbc_diffmode, y, Constant(cache.p))
     end
-    cache_collocation = __sparse_jacobian_cache(
-        Val(iip), jac_alg.nonbc_diffmode, sd_collocation,
-        loss_collocationₚ, resid_collocation, y)
 
-    J_bc = zero(init_jacobian(cache_bc))
-    J_c = zero(init_jacobian(cache_collocation))
     if J_full_band === nothing
         jac_prototype = vcat(J_bc, J_c)
     else
@@ -607,12 +669,12 @@ function __construct_nlproblem(
 
     jac = if iip
         @closure (J, u, p) -> __firk_mpoint_jacobian!(
-            J, J_c, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode, cache_bc,
-            cache_collocation, loss_bcₚ, loss_collocationₚ, resid_bc, resid_collocation, L)
+            J, J_c, u, bc_diffmode, nonbc_diffmode, cache_bc, cache_collocation,
+            loss_bc, loss_collocation, resid_bc, resid_collocation, L, cache.p)
     else
         @closure (u, p) -> __firk_mpoint_jacobian(
-            jac_prototype, J_c, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode,
-            cache_bc, cache_collocation, loss_bcₚ, loss_collocationₚ, L)
+            jac_prototype, J_c, u, bc_diffmode, nonbc_diffmode, cache_bc,
+            cache_collocation, loss_bc, loss_collocation, L, cache.p)
     end
 
     resid_prototype = vcat(resid_bc, resid_collocation)
@@ -628,31 +690,47 @@ function __construct_nlproblem(
     (; jac_alg) = cache.alg
     N = length(cache.mesh)
 
-    lossₚ = iip ? ((du, u) -> loss(du, u, cache.p)) : (u -> loss(u, cache.p))
-
     resid = vcat(@view(cache.bcresid_prototype[1:prod(cache.resid_size[1])]),
         __similar(y, cache.M * (N - 1)),
         @view(cache.bcresid_prototype[(prod(cache.resid_size[1]) + 1):end]))
     L = length(cache.bcresid_prototype)
 
-    sd = if jac_alg.diffmode isa AutoSparse
-        __sparsity_detection_alg(__generate_sparse_jacobian_prototype(
-            cache, cache.problem_type,
+    diffmode = if jac_alg.diffmode isa AutoSparse
+        colored_result = __generate_sparse_jacobian_prototype(cache, cache.problem_type,
             @view(cache.bcresid_prototype[1:prod(cache.resid_size[1])]),
             @view(cache.bcresid_prototype[(prod(cache.resid_size[1]) + 1):end]),
-            cache.M, N))
+            cache.M, N, jac_alg.diffmode)
+        partition = ifelse(
+            ADTypes.mode(jac_alg.bc_diffmode) isa ADTypes.ReverseMode, :row, :column)
+        constant_matrix_coloring = ifelse(
+            ADTypes.mode(jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode,
+            row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(jac_alg.diffmode);
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
+    else
+        jac_alg.diffmode
+    end
+
+    diffcache = if iip
+        DI.prepare_jacobian(loss, resid, diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss, diffmode, y, Constant(cache.p))
+    end
+
+    jac_prototype = if iip
+        DI.jacobian(loss, resid, diffcache, diffmode, y, Constant(cache.p))
     else
-        NoSparsityDetection()
+        DI.jacobian(loss, diffcache, diffmode, y, Constant(cache.p))
     end
-    diffcache = __sparse_jacobian_cache(Val(iip), jac_alg.diffmode, sd, lossₚ, resid, y)
-    jac_prototype = zero(init_jacobian(diffcache))
 
     jac = if iip
         @closure (J, u, p) -> __firk_2point_jacobian!(
-            J, u, jac_alg.diffmode, diffcache, lossₚ, resid)
+            J, u, jac_alg.diffmode, diffcache, loss, resid, cache.p)
     else
         @closure (u, p) -> __firk_2point_jacobian(
-            u, jac_prototype, jac_alg.diffmode, diffcache, lossₚ)
+            u, jac_prototype, jac_alg.diffmode, diffcache, loss, cache.p)
     end
 
     resid_prototype = copy(resid)
@@ -734,20 +812,21 @@ end
 
 function __firk_mpoint_jacobian!(
         J, _, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache, loss_bc::BC,
-        loss_collocation::C, resid_bc, resid_collocation, L::Int) where {BC, C}
-    sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, resid_bc, x)
-    sparse_jacobian!(@view(J[(L + 1):end, :]), nonbc_diffmode,
-        nonbc_diffcache, loss_collocation, resid_collocation, x)
+        loss_collocation::C, resid_bc, resid_collocation, L::Int, p) where {BC, C}
+    DI.jacobian!(
+        loss_bc, resid_bc, @view(J[1:L, :]), bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, resid_collocation, @view(J[(L + 1):end, :]),
+        nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     return nothing
 end
 
 function __firk_mpoint_jacobian!(J::AlmostBandedMatrix, J_c, x, bc_diffmode, nonbc_diffmode,
         bc_diffcache, nonbc_diffcache, loss_bc::BC, loss_collocation::C,
-        resid_bc, resid_collocation, L::Int) where {BC, C}
+        resid_bc, resid_collocation, L::Int, p) where {BC, C}
     J_bc = fillpart(J)
-    sparse_jacobian!(J_bc, bc_diffmode, bc_diffcache, loss_bc, resid_bc, x)
-    sparse_jacobian!(
-        J_c, nonbc_diffmode, nonbc_diffcache, loss_collocation, resid_collocation, x)
+    DI.jacobian!(loss_bc, resid_bc, J_bc, bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, resid_collocation, J_c,
+        nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     exclusive_bandpart(J) .= J_c
     finish_part_setindex!(J)
     return nothing
@@ -755,30 +834,30 @@ end
 
 function __firk_mpoint_jacobian(
         J, _, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache,
-        loss_bc::BC, loss_collocation::C, L::Int) where {BC, C}
-    sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, x)
-    sparse_jacobian!(
-        @view(J[(L + 1):end, :]), nonbc_diffmode, nonbc_diffcache, loss_collocation, x)
+        loss_bc::BC, loss_collocation::C, L::Int, p) where {BC, C}
+    DI.jacobian!(loss_bc, @view(J[1:L, :]), bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, @view(J[(L + 1):end, :]),
+        nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     return J
 end
 
 function __firk_mpoint_jacobian(
         J::AlmostBandedMatrix, J_c, x, bc_diffmode, nonbc_diffmode, bc_diffcache,
-        nonbc_diffcache, loss_bc::BC, loss_collocation::C, L::Int) where {BC, C}
+        nonbc_diffcache, loss_bc::BC, loss_collocation::C, L::Int, p) where {BC, C}
     J_bc = fillpart(J)
-    sparse_jacobian!(J_bc, bc_diffmode, bc_diffcache, loss_bc, x)
-    sparse_jacobian!(J_c, nonbc_diffmode, nonbc_diffcache, loss_collocation, x)
+    DI.jacobian!(loss_bc, J_bc, bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, J_c, nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     exclusive_bandpart(J) .= J_c
     finish_part_setindex!(J)
     return J
 end
 
-function __firk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid) where {L}
-    sparse_jacobian!(J, diffmode, diffcache, loss_fn, resid, x)
+function __firk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid, p) where {L}
+    DI.jacobian!(loss_fn, resid, J, diffcache, diffmode, x, Constant(p))
     return J
 end
 
-function __firk_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L) where {L}
-    sparse_jacobian!(J, diffmode, diffcache, loss_fn, x)
+function __firk_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L, p) where {L}
+    DI.jacobian!(loss_fn, J, diffcache, diffmode, x, Constant(p))
     return J
 end
diff --git a/lib/BoundaryValueDiffEqFIRK/src/sparse_jacobians.jl b/lib/BoundaryValueDiffEqFIRK/src/sparse_jacobians.jl
index 9d1c5349..de8137b3 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/sparse_jacobians.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/sparse_jacobians.jl
@@ -12,27 +12,34 @@ If the problem is a TwoPointBVProblem, then this is the complete Jacobian, else
 computes the sparse part excluding the contributions from the boundary conditions.
 """
 function __generate_sparse_jacobian_prototype(
-        ::FIRKCacheNested, ::StandardBVProblem, ya, yb, M, N)
+        ::FIRKCacheNested, ::StandardBVProblem, ya, yb, M, N, ad)
     fast_scalar_indexing(ya) ||
         error("Sparse Jacobians are only supported for Fast Scalar Index-able Arrays")
     J_c = BandedMatrix(Ones{eltype(ya)}(M * (N - 1), M * N), (1, 2M - 1))
-    return ColoredMatrix(J_c, matrix_colors(J_c'), matrix_colors(J_c))
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
+    result = coloring(J_c, problem, algo)
+    return result
 end
 
 function __generate_sparse_jacobian_prototype(
-        ::FIRKCacheNested, ::TwoPointBVProblem, ya, yb, M, N)
+        ::FIRKCacheNested, ::TwoPointBVProblem, ya, yb, M, N, ad)
     fast_scalar_indexing(ya) ||
         error("Sparse Jacobians are only supported for Fast Scalar Index-able Arrays")
     J₁ = length(ya) + length(yb) + M * (N - 1)
     J₂ = M * N
     J = BandedMatrix(Ones{eltype(ya)}(J₁, J₂), (M + 1, M + 1))
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
     # for underdetermined systems we don't have banded qr implemented. use sparse
-    J₁ < J₂ && return ColoredMatrix(sparse(J), matrix_colors(J'), matrix_colors(J))
-    return ColoredMatrix(J, matrix_colors(J'), matrix_colors(J))
+    J₁ < J₂ && return coloring(sparse(J), problem, algo)
+    return coloring(J, problem, algo)
 end
 
 function __generate_sparse_jacobian_prototype(
-        cache::FIRKCacheExpand, ::StandardBVProblem, ya, yb, M, N)
+        cache::FIRKCacheExpand, ::StandardBVProblem, ya, yb, M, N, ad)
     (; stage) = cache
 
     # Get number of nonzeros
@@ -65,20 +72,14 @@ function __generate_sparse_jacobian_prototype(
     # Create sparse matrix from Is and Js
     J_c = _sparse_like(Is, Js, ya, row_size, row_size + M)
 
-    col_colorvec = Vector{Int}(undef, size(J_c, 2))
-    for i in eachindex(col_colorvec)
-        col_colorvec[i] = mod1(i, (2 * M * (stage + 1)) + M)
-    end
-    row_colorvec = Vector{Int}(undef, size(J_c, 1))
-    for i in eachindex(row_colorvec)
-        row_colorvec[i] = mod1(i, (2 * M * (stage + 1)) + M)
-    end
-
-    return ColoredMatrix(J_c, row_colorvec, col_colorvec)
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
+    return coloring(J_c, problem, algo)
 end
 
 function __generate_sparse_jacobian_prototype(
-        cache::FIRKCacheExpand, ::TwoPointBVProblem, ya, yb, M, N)
+        cache::FIRKCacheExpand, ::TwoPointBVProblem, ya, yb, M, N, ad)
     (; stage) = cache
 
     # Get number of nonzeros
@@ -129,5 +130,9 @@ function __generate_sparse_jacobian_prototype(
 
     # Create sparse matrix from Is and Js
     J = _sparse_like(Is, Js, ya, row_size + length(ya) + length(yb), row_size + M)
-    return ColoredMatrix(J, matrix_colors(J'), matrix_colors(J))
+
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
+    return coloring(J, problem, algo)
 end
diff --git a/lib/BoundaryValueDiffEqFIRK/src/utils.jl b/lib/BoundaryValueDiffEqFIRK/src/utils.jl
index 371516f1..d3b5a274 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/utils.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/utils.jl
@@ -1,4 +1,5 @@
-function __append_similar!(x::AbstractVector{<:AbstractArray}, n, _, TU::FIRKTableau{false})
+function BoundaryValueDiffEqCore.__append_similar!(
+        x::AbstractVector{<:AbstractArray}, n, _, TU::FIRKTableau{false})
     (; s) = TU
     N = (n - 1) * (s + 1) + 1 - length(x)
     N == 0 && return x
@@ -7,7 +8,7 @@ function __append_similar!(x::AbstractVector{<:AbstractArray}, n, _, TU::FIRKTab
     return x
 end
 
-function __append_similar!(
+function BoundaryValueDiffEqCore.__append_similar!(
         x::AbstractVector{<:MaybeDiffCache}, n, M, TU::FIRKTableau{false})
     (; s) = TU
     N = (n - 1) * (s + 1) + 1 - length(x)
@@ -19,7 +20,8 @@ function __append_similar!(
     return x
 end
 
-function __append_similar!(x::AbstractVectorOfArray, n, M, TU::FIRKTableau{false})
+function BoundaryValueDiffEqCore.__append_similar!(
+        x::AbstractVectorOfArray, n, M, TU::FIRKTableau{false})
     (; s) = TU
     N = (n - 1) * (s + 1) + 1 - length(x)
     N == 0 && return x
@@ -28,7 +30,8 @@ function __append_similar!(x::AbstractVectorOfArray, n, M, TU::FIRKTableau{false
     return x
 end
 
-function __append_similar!(x::AbstractVectorOfArray, n, M, TU::FIRKTableau{true})
+function BoundaryValueDiffEqCore.__append_similar!(
+        x::AbstractVectorOfArray, n, M, TU::FIRKTableau{true})
     (; s) = TU
     N = n - length(x)
     N == 0 && return x
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index 5876991c..66975fcf 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -180,8 +180,10 @@ end
         @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5)
             @time sim = test_convergence(
                 dts, prob, lobattoIIIa_solver(Val(stage)); abstol = 1e-8)
-            if (stage == 5) || (((i == 9) || (i == 10)) && stage == 4)
-                @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
+            if (stage == 5)
+                @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=0.4
+            elseif (stage == 4)
+                @test sim.𝒪est[:final]≈2 * stage - 2 atol=0.5
             else
                 @test sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             end
@@ -193,7 +195,7 @@ end
             if (stage == 5) || (stage == 4 && i == 10)
                 @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             elseif stage == 4
-                @test sim.𝒪est[:final]≈2 * stage - 2 atol=0.5
+                @test sim.𝒪est[:final]≈2 * stage - 2 atol=0.7
             else
                 @test sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             end
@@ -202,7 +204,7 @@ end
         @testset "LobattoIIIc$stage" for stage in (2, 3, 4, 5)
             @time sim = test_convergence(
                 dts, prob, lobattoIIIc_solver(Val(stage)); abstol = 1e-8, reltol = 1e-8)
-            if (i == 3 && stage == 4) || (i == 4 && stage == 4)
+            if stage == 4
                 @test sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             elseif first(sim.errors[:final]) < 1e-12
                 @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
diff --git a/lib/BoundaryValueDiffEqMIRK/Project.toml b/lib/BoundaryValueDiffEqMIRK/Project.toml
index 4287aa6a..a0ddfdbd 100644
--- a/lib/BoundaryValueDiffEqMIRK/Project.toml
+++ b/lib/BoundaryValueDiffEqMIRK/Project.toml
@@ -1,6 +1,6 @@
 name = "BoundaryValueDiffEqMIRK"
 uuid = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
-version = "1.2.0"
+version = "1.3.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -10,12 +10,12 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
@@ -24,7 +24,8 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 
 [compat]
 ADTypes = "1.9"
@@ -36,6 +37,7 @@ BoundaryValueDiffEqCore = "1.1.0"
 ConcreteStructs = "0.2.3"
 DiffEqBase = "6.158.3"
 DiffEqDevTools = "2.44"
+DifferentiationInterface = "0.6.24"
 FastAlmostBandedMatrices = "0.1.4"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.38"
@@ -44,7 +46,6 @@ InteractiveUtils = "<0.0.1, 1"
 JET = "0.9.12"
 LinearAlgebra = "1.10"
 LinearSolve = "2.36.2"
-Logging = "1.10"
 OrdinaryDiffEq = "6.90.1"
 PreallocationTools = "0.4.24"
 PrecompileTools = "1.2"
@@ -56,7 +57,8 @@ Reexport = "1.2"
 SciMLBase = "2.60.0"
 Setfield = "1.1.1"
 SparseArrays = "1.10"
-SparseDiffTools = "2.23"
+SparseConnectivityTracer = "0.6.9"
+SparseMatrixColorings = "0.4.10"
 StaticArrays = "1.9.8"
 Test = "1.10"
 julia = "1.10"
diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
index f271c207..dba5c420 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
@@ -1,44 +1,51 @@
 module BoundaryValueDiffEqMIRK
 
-import PrecompileTools: @compile_workload, @setup_workload
-
-using ADTypes, Adapt, ArrayInterface, BoundaryValueDiffEqCore, DiffEqBase, ForwardDiff,
-      LinearAlgebra, Preferences, RecursiveArrayTools, Reexport, SciMLBase, Setfield,
-      SparseDiffTools
-
+using ADTypes
+using ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
+using BandedMatrices: BandedMatrix, Ones
+using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
+                               recursive_flatten, recursive_flatten!, recursive_unflatten!,
+                               __concrete_nonlinearsolve_algorithm, diff!,
+                               __FastShortcutBVPCompatibleNonlinearPolyalg,
+                               __FastShortcutBVPCompatibleNLLSPolyalg, EvalSol,
+                               concrete_jacobian_algorithm, eval_bc_residual,
+                               eval_bc_residual!, get_tmp, __maybe_matmul!,
+                               __append_similar!, __extract_problem_details,
+                               __initial_guess, __maybe_allocate_diffcache,
+                               __restructure_sol, __get_bcresid_prototype, __similar, __vec,
+                               __vec_f, __vec_f!, __vec_bc, __vec_bc!,
+                               recursive_flatten_twopoint!, __internal_nlsolve_problem,
+                               __extract_mesh, __extract_u0, __has_initial_guess,
+                               __initial_guess_length, __initial_guess_on_mesh,
+                               __flatten_initial_guess, __build_solution, __Fix3,
+                               _sparse_like, get_dense_ad, _sparse_like
+
+using ConcreteStructs: @concrete
+using DiffEqBase: DiffEqBase
+using DifferentiationInterface: DifferentiationInterface, Constant, prepare_jacobian
+using FastAlmostBandedMatrices: AlmostBandedMatrix, fillpart, exclusive_bandpart,
+                                finish_part_setindex!
+using FastClosures: @closure
+using ForwardDiff: ForwardDiff, pickchunksize, Dual
+using LinearAlgebra
+using RecursiveArrayTools: AbstractVectorOfArray, DiffEqArray, VectorOfArray, recursivecopy,
+                           recursivefill!
+using SciMLBase: SciMLBase, AbstractDiffEqInterpolation, StandardBVProblem, __solve,
+                 _unwrap_val
+using Setfield: @set!
+using Reexport: @reexport
 using PreallocationTools: PreallocationTools, DiffCache
+using PrecompileTools: @compile_workload, @setup_workload
+using Preferences: Preferences
+using SparseArrays: sparse
+using SparseConnectivityTracer: SparseConnectivityTracer
+using SparseMatrixColorings: ColoringProblem, GreedyColoringAlgorithm, sparsity_pattern,
+                             ConstantColoringAlgorithm, row_colors, column_colors, coloring,
+                             LargestFirst
+
+const DI = DifferentiationInterface
 
-# Special Matrix Types
-using BandedMatrices, FastAlmostBandedMatrices, SparseArrays
-
-import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
-                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                                __concrete_nonlinearsolve_algorithm, diff!,
-                                __FastShortcutBVPCompatibleNonlinearPolyalg,
-                                __FastShortcutBVPCompatibleNLLSPolyalg, __restructure_sol,
-                                concrete_jacobian_algorithm, eval_bc_residual,
-                                eval_bc_residual!, get_tmp, __maybe_matmul!,
-                                __append_similar!, __extract_problem_details,
-                                __initial_guess, __maybe_allocate_diffcache,
-                                __get_bcresid_prototype, __similar, __vec, __vec_f,
-                                __vec_f!, __vec_bc, __vec_bc!, recursive_flatten_twopoint!,
-                                __internal_nlsolve_problem, __extract_mesh, __extract_u0,
-                                __has_initial_guess, __initial_guess_length, EvalSol,
-                                __initial_guess_on_mesh, __flatten_initial_guess,
-                                __build_solution, __Fix3, __sparse_jacobian_cache,
-                                __sparsity_detection_alg, _sparse_like, ColoredMatrix
-
-import ADTypes: AbstractADType
-import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
-import ConcreteStructs: @concrete
-import DiffEqBase: solve
-import FastClosures: @closure
-import ForwardDiff: ForwardDiff, pickchunksize, Dual
-import Logging
-import RecursiveArrayTools: ArrayPartition, DiffEqArray
-import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
-
-@reexport using ADTypes, BoundaryValueDiffEqCore, SparseDiffTools, SciMLBase
+@reexport using ADTypes, BoundaryValueDiffEqCore, SciMLBase
 
 include("types.jl")
 include("algorithms.jl")
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index b69a8450..53a36683 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -319,39 +319,64 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
     L = length(resid_bc)
     resid_collocation = __similar(y, cache.M * (N - 1))
 
-    loss_bcₚ = (iip ? __Fix3 : Base.Fix2)(loss_bc, cache.p)
-    loss_collocationₚ = (iip ? __Fix3 : Base.Fix2)(loss_collocation, cache.p)
+    bc_diffmode = if jac_alg.bc_diffmode isa AutoSparse
+        AutoSparse(jac_alg.bc_diffmode;
+            sparsity_detector = SparseConnectivityTracer.TracerSparsityDetector(),
+            coloring_algorithm = GreedyColoringAlgorithm())
+    else
+        jac_alg.bc_diffmode
+    end
 
-    sd_bc = jac_alg.bc_diffmode isa AutoSparse ? SymbolicsSparsityDetection() :
-            NoSparsityDetection()
-    cache_bc = __sparse_jacobian_cache(
-        Val(iip), jac_alg.bc_diffmode, sd_bc, loss_bcₚ, resid_bc, y)
+    cache_bc = if iip
+        DI.prepare_jacobian(loss_bc, resid_bc, bc_diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss_bc, bc_diffmode, y, Constant(cache.p))
+    end
 
-    sd_collocation = if jac_alg.nonbc_diffmode isa AutoSparse
+    nonbc_diffmode = if jac_alg.nonbc_diffmode isa AutoSparse
         if L < cache.M
             # For underdetermined problems we use sparse since we don't have banded qr
-            colored_matrix = __generate_sparse_jacobian_prototype(
-                cache, cache.problem_type, y, y, cache.M, N)
             J_full_band = nothing
-            __sparsity_detection_alg(ColoredMatrix(
-                sparse(colored_matrix.M), colored_matrix.row_colorvec,
-                colored_matrix.col_colorvec))
+            colored_result = __generate_sparse_jacobian_prototype(
+                cache, cache.problem_type, y, y, cache.M, N, jac_alg.nonbc_diffmode)
         else
             J_full_band = BandedMatrix(Ones{eltype(y)}(L + cache.M * (N - 1), cache.M * N),
                 (L + 1, cache.M + max(cache.M - L, 0)))
-            __sparsity_detection_alg(__generate_sparse_jacobian_prototype(
-                cache, cache.problem_type, y, y, cache.M, N))
+            colored_result = __generate_sparse_jacobian_prototype(
+                cache, cache.problem_type, y, y, cache.M, N, jac_alg.nonbc_diffmode)
         end
+        is_diffmode_reverse = ADTypes.mode(jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode
+        partition = ifelse(is_diffmode_reverse, :row, :column)
+        constant_matrix_coloring = ifelse(is_diffmode_reverse, row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(jac_alg.nonbc_diffmode);
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
     else
         J_full_band = nothing
-        NoSparsityDetection()
+        jac_alg.nonbc_diffmode
+    end
+
+    cache_collocation = if iip
+        DI.prepare_jacobian(
+            loss_collocation, resid_collocation, nonbc_diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss_collocation, nonbc_diffmode, y, Constant(cache.p))
+    end
+
+    J_bc = if iip
+        DI.jacobian(loss_bc, resid_bc, cache_bc, bc_diffmode, y, Constant(cache.p))
+    else
+        DI.jacobian(loss_bc, cache_bc, bc_diffmode, y, Constant(cache.p))
+    end
+    J_c = if iip
+        DI.jacobian(loss_collocation, resid_collocation, cache_collocation,
+            nonbc_diffmode, y, Constant(cache.p))
+    else
+        DI.jacobian(
+            loss_collocation, cache_collocation, nonbc_diffmode, y, Constant(cache.p))
     end
-    cache_collocation = __sparse_jacobian_cache(
-        Val(iip), jac_alg.nonbc_diffmode, sd_collocation,
-        loss_collocationₚ, resid_collocation, y)
 
-    J_bc = zero(init_jacobian(cache_bc))
-    J_c = zero(init_jacobian(cache_collocation))
     if J_full_band === nothing
         jac_prototype = vcat(J_bc, J_c)
     else
@@ -360,12 +385,12 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
 
     jac = if iip
         @closure (J, u, p) -> __mirk_mpoint_jacobian!(
-            J, J_c, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode, cache_bc,
-            cache_collocation, loss_bcₚ, loss_collocationₚ, resid_bc, resid_collocation, L)
+            J, J_c, u, bc_diffmode, nonbc_diffmode, cache_bc, cache_collocation,
+            loss_bc, loss_collocation, resid_bc, resid_collocation, L, cache.p)
     else
         @closure (u, p) -> __mirk_mpoint_jacobian(
-            jac_prototype, J_c, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode,
-            cache_bc, cache_collocation, loss_bcₚ, loss_collocationₚ, L)
+            jac_prototype, J_c, u, bc_diffmode, nonbc_diffmode, cache_bc,
+            cache_collocation, loss_bc, loss_collocation, L, cache.p)
     end
 
     resid_prototype = vcat(resid_bc, resid_collocation)
@@ -377,20 +402,21 @@ end
 
 function __mirk_mpoint_jacobian!(
         J, _, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache, loss_bc::BC,
-        loss_collocation::C, resid_bc, resid_collocation, L::Int) where {BC, C}
-    sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, resid_bc, x)
-    sparse_jacobian!(@view(J[(L + 1):end, :]), nonbc_diffmode,
-        nonbc_diffcache, loss_collocation, resid_collocation, x)
+        loss_collocation::C, resid_bc, resid_collocation, L::Int, p) where {BC, C}
+    DI.jacobian!(
+        loss_bc, resid_bc, @view(J[1:L, :]), bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, resid_collocation, @view(J[(L + 1):end, :]),
+        nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     return nothing
 end
 
 function __mirk_mpoint_jacobian!(J::AlmostBandedMatrix, J_c, x, bc_diffmode, nonbc_diffmode,
         bc_diffcache, nonbc_diffcache, loss_bc::BC, loss_collocation::C,
-        resid_bc, resid_collocation, L::Int) where {BC, C}
+        resid_bc, resid_collocation, L::Int, p) where {BC, C}
     J_bc = fillpart(J)
-    sparse_jacobian!(J_bc, bc_diffmode, bc_diffcache, loss_bc, resid_bc, x)
-    sparse_jacobian!(
-        J_c, nonbc_diffmode, nonbc_diffcache, loss_collocation, resid_collocation, x)
+    DI.jacobian!(loss_bc, resid_bc, J_bc, bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, resid_collocation, J_c,
+        nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     exclusive_bandpart(J) .= J_c
     finish_part_setindex!(J)
     return nothing
@@ -398,19 +424,19 @@ end
 
 function __mirk_mpoint_jacobian(
         J, _, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache,
-        loss_bc::BC, loss_collocation::C, L::Int) where {BC, C}
-    sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, x)
-    sparse_jacobian!(
-        @view(J[(L + 1):end, :]), nonbc_diffmode, nonbc_diffcache, loss_collocation, x)
+        loss_bc::BC, loss_collocation::C, L::Int, p) where {BC, C}
+    DI.jacobian!(loss_bc, @view(J[1:L, :]), bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, @view(J[(L + 1):end, :]),
+        nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     return J
 end
 
 function __mirk_mpoint_jacobian(
         J::AlmostBandedMatrix, J_c, x, bc_diffmode, nonbc_diffmode, bc_diffcache,
-        nonbc_diffcache, loss_bc::BC, loss_collocation::C, L::Int) where {BC, C}
+        nonbc_diffcache, loss_bc::BC, loss_collocation::C, L::Int, p) where {BC, C}
     J_bc = fillpart(J)
-    sparse_jacobian!(J_bc, bc_diffmode, bc_diffcache, loss_bc, x)
-    sparse_jacobian!(J_c, nonbc_diffmode, nonbc_diffcache, loss_collocation, x)
+    DI.jacobian!(loss_bc, J_bc, bc_diffcache, bc_diffmode, x, Constant(p))
+    DI.jacobian!(loss_collocation, J_c, nonbc_diffcache, nonbc_diffmode, x, Constant(p))
     exclusive_bandpart(J) .= J_c
     finish_part_setindex!(J)
     return J
@@ -421,31 +447,47 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
     (; nlsolve, jac_alg) = cache.alg
     N = length(cache.mesh)
 
-    lossₚ = iip ? ((du, u) -> loss(du, u, cache.p)) : (u -> loss(u, cache.p))
-
     resid = vcat(@view(cache.bcresid_prototype[1:prod(cache.resid_size[1])]),
         __similar(y, cache.M * (N - 1)),
         @view(cache.bcresid_prototype[(prod(cache.resid_size[1]) + 1):end]))
     L = length(cache.bcresid_prototype)
 
-    sd = if jac_alg.diffmode isa AutoSparse
-        __sparsity_detection_alg(__generate_sparse_jacobian_prototype(
-            cache, cache.problem_type,
+    diffmode = if jac_alg.diffmode isa AutoSparse
+        colored_result = __generate_sparse_jacobian_prototype(cache, cache.problem_type,
             @view(cache.bcresid_prototype[1:prod(cache.resid_size[1])]),
             @view(cache.bcresid_prototype[(prod(cache.resid_size[1]) + 1):end]),
-            cache.M, N))
+            cache.M, N, jac_alg.diffmode)
+        partition = ifelse(
+            ADTypes.mode(jac_alg.bc_diffmode) isa ADTypes.ReverseMode, :row, :column)
+        constant_matrix_coloring = ifelse(
+            ADTypes.mode(jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode,
+            row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(jac_alg.diffmode);
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
+    else
+        jac_alg.diffmode
+    end
+
+    diffcache = if iip
+        DI.prepare_jacobian(loss, resid, diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss, diffmode, y, Constant(cache.p))
+    end
+
+    jac_prototype = if iip
+        DI.jacobian(loss, resid, diffcache, diffmode, y, Constant(cache.p)) #zero(init_jacobian(diffcache))
     else
-        NoSparsityDetection()
+        DI.jacobian(loss, diffcache, diffmode, y, Constant(cache.p))
     end
-    diffcache = __sparse_jacobian_cache(Val(iip), jac_alg.diffmode, sd, lossₚ, resid, y)
-    jac_prototype = zero(init_jacobian(diffcache))
 
     jac = if iip
         @closure (J, u, p) -> __mirk_2point_jacobian!(
-            J, u, jac_alg.diffmode, diffcache, lossₚ, resid)
+            J, u, jac_alg.diffmode, diffcache, loss, resid, p)
     else
         @closure (u, p) -> __mirk_2point_jacobian(
-            u, jac_prototype, jac_alg.diffmode, diffcache, lossₚ)
+            u, jac_prototype, jac_alg.diffmode, diffcache, loss, p)
     end
 
     resid_prototype = copy(resid)
@@ -454,12 +496,12 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
     return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
 end
 
-function __mirk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid) where {L}
-    sparse_jacobian!(J, diffmode, diffcache, loss_fn, resid, x)
+function __mirk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid, p) where {L}
+    DI.jacobian!(loss_fn, resid, J, diffcache, diffmode, x, Constant(p))
     return J
 end
 
-function __mirk_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L) where {L}
-    sparse_jacobian!(J, diffmode, diffcache, loss_fn, x)
+function __mirk_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L, p) where {L}
+    DI.jacobian!(loss_fn, J, diffcache, diffmode, x, Constant(p))
     return J
 end
diff --git a/lib/BoundaryValueDiffEqMIRK/src/sparse_jacobians.jl b/lib/BoundaryValueDiffEqMIRK/src/sparse_jacobians.jl
index 4b28ebcb..051e5201 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/sparse_jacobians.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/sparse_jacobians.jl
@@ -10,26 +10,33 @@ coloring.
 If the problem is a TwoPointBVProblem, then this is the complete Jacobian, else it only
 computes the sparse part excluding the contributions from the boundary conditions.
 """
-function __generate_sparse_jacobian_prototype(cache::MIRKCache, ya, yb, M, N)
-    return __generate_sparse_jacobian_prototype(cache, cache.problem_type, ya, yb, M, N)
+function __generate_sparse_jacobian_prototype(cache::MIRKCache, ya, yb, M, N, ad)
+    return __generate_sparse_jacobian_prototype(cache, cache.problem_type, ya, yb, M, N, ad)
 end
 
 function __generate_sparse_jacobian_prototype(
-        ::MIRKCache, ::StandardBVProblem, ya, yb, M, N)
+        ::MIRKCache, ::StandardBVProblem, ya, yb, M, N, ad)
     fast_scalar_indexing(ya) ||
         error("Sparse Jacobians are only supported for Fast Scalar Index-able Arrays")
     J_c = BandedMatrix(Ones{eltype(ya)}(M * (N - 1), M * N), (1, 2M - 1))
-    return ColoredMatrix(J_c, matrix_colors(J_c'), matrix_colors(J_c))
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
+    result = coloring(J_c, problem, algo)
+    return result
 end
 
 function __generate_sparse_jacobian_prototype(
-        ::MIRKCache, ::TwoPointBVProblem, ya, yb, M, N)
+        ::MIRKCache, ::TwoPointBVProblem, ya, yb, M, N, ad)
     fast_scalar_indexing(ya) ||
         error("Sparse Jacobians are only supported for Fast Scalar Index-able Arrays")
     J₁ = length(ya) + length(yb) + M * (N - 1)
     J₂ = M * N
     J = BandedMatrix(Ones{eltype(ya)}(J₁, J₂), (M + 1, M + 1))
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
     # for underdetermined systems we don't have banded qr implemented. use sparse
-    J₁ < J₂ && return ColoredMatrix(sparse(J), matrix_colors(J'), matrix_colors(J))
-    return ColoredMatrix(J, matrix_colors(J'), matrix_colors(J))
+    J₁ < J₂ && return coloring(sparse(J), problem, algo)
+    return coloring(J, problem, algo)
 end
diff --git a/lib/BoundaryValueDiffEqMIRKN/Project.toml b/lib/BoundaryValueDiffEqMIRKN/Project.toml
index 9e7f0a5d..bf96f4e1 100644
--- a/lib/BoundaryValueDiffEqMIRKN/Project.toml
+++ b/lib/BoundaryValueDiffEqMIRKN/Project.toml
@@ -1,7 +1,7 @@
 name = "BoundaryValueDiffEqMIRKN"
 uuid = "9255f1d6-53bf-473e-b6bd-23f1ff009da4"
 authors = ["Qingyu Qu <erikqqy123@gmail.com>"]
-version = "1.0.0"
+version = "1.1.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -11,13 +11,13 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
@@ -26,7 +26,8 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 
 [compat]
 ADTypes = "1.9"
@@ -38,6 +39,7 @@ BoundaryValueDiffEqCore = "1.0.0"
 ConcreteStructs = "0.2.3"
 DiffEqBase = "6.158.3"
 DiffEqDevTools = "2.44"
+DifferentiationInterface = "0.6.24"
 FastAlmostBandedMatrices = "0.1.4"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.38"
@@ -46,8 +48,7 @@ InteractiveUtils = "<0.0.1, 1"
 JET = "0.9.12"
 LineSearch = "0.1.3"
 LinearAlgebra = "1.10"
-LinearSolve = "2.21"
-Logging = "1.10"
+LinearSolve = "2.37.0"
 OrdinaryDiffEq = "6.90.1"
 PreallocationTools = "0.4.24"
 PrecompileTools = "1.2"
@@ -59,7 +60,8 @@ Reexport = "1.2"
 SciMLBase = "2.60.0"
 Setfield = "1.1.1"
 SparseArrays = "1.10"
-SparseDiffTools = "2.23"
+SparseConnectivityTracer = "0.6.9"
+SparseMatrixColorings = "0.4.10"
 StaticArrays = "1.9.8"
 Test = "1.10"
 julia = "1.10"
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
index aec4d55f..3ad2605a 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
@@ -1,46 +1,51 @@
 module BoundaryValueDiffEqMIRKN
 
-import PrecompileTools: @compile_workload, @setup_workload
-
-using ADTypes, Adapt, ArrayInterface, BoundaryValueDiffEqCore, DiffEqBase, ForwardDiff,
-      LinearAlgebra, Preferences, RecursiveArrayTools, Reexport, SciMLBase, Setfield,
-      SparseDiffTools
-
+using ADTypes: ADTypes, AutoSparse, AutoForwardDiff
+using ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
+using BandedMatrices: BandedMatrix, Ones
+using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
+                               recursive_flatten, recursive_flatten!, recursive_unflatten!,
+                               __concrete_nonlinearsolve_algorithm, diff!, EvalSol,
+                               __FastShortcutBVPCompatibleNonlinearPolyalg,
+                               __FastShortcutBVPCompatibleNLLSPolyalg, eval_bc_residual,
+                               eval_bc_residual!, get_tmp, __maybe_matmul!,
+                               __append_similar!, __extract_problem_details,
+                               __initial_guess, __maybe_allocate_diffcache,
+                               __restructure_sol, __get_bcresid_prototype, __similar, __vec,
+                               __vec_f, __vec_f!, __vec_bc, __vec_bc!, __vec_so_bc!,
+                               __vec_so_bc, recursive_flatten_twopoint!,
+                               __internal_nlsolve_problem, __extract_mesh, __extract_u0,
+                               __has_initial_guess, __initial_guess_length,
+                               __initial_guess_on_mesh, __flatten_initial_guess,
+                               __build_solution, __Fix3, _sparse_like, __default_sparse_ad,
+                               __default_nonsparse_ad, get_dense_ad
+
+using ConcreteStructs: @concrete
+using DiffEqBase: DiffEqBase
+using DifferentiationInterface: DifferentiationInterface, Constant
+using FastAlmostBandedMatrices: AlmostBandedMatrix, fillpart, exclusive_bandpart,
+                                finish_part_setindex!
+using FastClosures: @closure
+using ForwardDiff: ForwardDiff, pickchunksize, Dual
+using LinearAlgebra
 using PreallocationTools: PreallocationTools, DiffCache
-
-# Special Matrix Types
-using BandedMatrices, FastAlmostBandedMatrices, SparseArrays
-
-import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
-                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                                __concrete_nonlinearsolve_algorithm, diff!,
-                                __FastShortcutBVPCompatibleNonlinearPolyalg,
-                                __FastShortcutBVPCompatibleNLLSPolyalg, eval_bc_residual,
-                                eval_bc_residual!, get_tmp, __maybe_matmul!, EvalSol,
-                                __append_similar!, __extract_problem_details,
-                                __initial_guess, __maybe_allocate_diffcache,
-                                __restructure_sol, __get_bcresid_prototype, __similar,
-                                __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
-                                recursive_flatten_twopoint!, __internal_nlsolve_problem,
-                                __extract_mesh, __extract_u0, __has_initial_guess,
-                                __initial_guess_length, __initial_guess_on_mesh,
-                                __flatten_initial_guess, __build_solution, __Fix3,
-                                __sparse_jacobian_cache, __sparsity_detection_alg,
-                                _sparse_like, ColoredMatrix, __default_sparse_ad,
-                                __default_nonsparse_ad
-
-import ADTypes: AbstractADType
-import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
-import ConcreteStructs: @concrete
-import DiffEqBase: solve
-import FastClosures: @closure
-import ForwardDiff: ForwardDiff, pickchunksize
-import Logging
-import RecursiveArrayTools: ArrayPartition, DiffEqArray
-import SciMLBase: AbstractDiffEqInterpolation, AbstractBVProblem,
-                  StandardSecondOrderBVProblem, StandardBVProblem, __solve, _unwrap_val
-
-@reexport using ADTypes, BoundaryValueDiffEqCore, SparseDiffTools, SciMLBase
+using PrecompileTools: @compile_workload, @setup_workload
+using Preferences: Preferences
+using RecursiveArrayTools: AbstractVectorOfArray, VectorOfArray, recursivecopy,
+                           ArrayPartition
+using Reexport: @reexport
+using SciMLBase: SciMLBase, AbstractDiffEqInterpolation, AbstractBVProblem,
+                 StandardSecondOrderBVProblem, StandardBVProblem, __solve, _unwrap_val
+using Setfield: @set!, @set
+using SparseArrays: sparse
+using SparseConnectivityTracer: SparseConnectivityTracer
+using SparseMatrixColorings: ColoringProblem, GreedyColoringAlgorithm,
+                             ConstantColoringAlgorithm, row_colors, column_colors, coloring,
+                             LargestFirst
+
+const DI = DifferentiationInterface
+
+@reexport using ADTypes, BoundaryValueDiffEqCore, SciMLBase
 
 include("utils.jl")
 include("types.jl")
@@ -51,12 +56,6 @@ include("collocation.jl")
 include("mirkn_tableaus.jl")
 include("interpolation.jl")
 
-function __solve(
-        prob::AbstractBVProblem, alg::BoundaryValueDiffEqAlgorithm, args...; kwargs...)
-    cache = init(prob, alg, args...; kwargs...)
-    return solve!(cache)
-end
-
 export MIRKN4, MIRKN6
 export BVPJacobianAlgorithm
 
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
index 3909b0f2..114a2fd3 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
@@ -39,9 +39,9 @@ function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN;
     # Don't flatten this here, since we need to expand it later if needed
     y₀ = __initial_guess_on_mesh(prob, prob.u0, Nig, prob.p, false)
     chunksize = pickchunksize(M * (2 * Nig - 2))
-    __alloc = @closure x -> __maybe_allocate_diffcache(vec(x), chunksize, alg.jac_alg)
+    __alloc = @closure x -> __maybe_allocate_diffcache(vec(zero(x)), chunksize, alg.jac_alg)
 
-    y = __alloc.(copy.(y₀))
+    y = __alloc.(copy.(y₀.u))
     fᵢ_cache = __alloc(similar(X))
     fᵢ₂_cache = __alloc(similar(X))
     stage = alg_stage(alg)
@@ -50,7 +50,7 @@ function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN;
                   for _ in 1:Nig]
 
     residual = if iip
-        __alloc.(copy.(@view(y₀[1:end])))
+        __alloc.(copy.(@view(y₀.u[1:end])))
     else
         nothing
     end
@@ -96,14 +96,16 @@ function SciMLBase.solve!(cache::MIRKNCache{iip, T}) where {iip, T}
     nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve)
     sol_nlprob = __solve(nlprob, nlsolve_alg; kwargs..., alias_u0 = true)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
-    solu = ArrayPartition.(cache.y₀[1:length(mesh)], cache.y₀[(length(mesh) + 1):end])
+    L = length(mesh)
+    solu = ArrayPartition.(cache.y₀.u[1:L], cache.y₀.u[(L + 1):end])
     return SciMLBase.build_solution(
         prob, cache.alg, mesh, solu; retcode = sol_nlprob.retcode)
 end
 
-function __construct_nlproblem(cache::MIRKNCache{iip}, y::AbstractVector) where {iip}
+function __construct_nlproblem(cache::MIRKNCache{iip}, y) where {iip}
     (; alg) = cache
     pt = cache.problem_type
+
     loss = if iip
         @closure (du, u, p) -> __mirkn_loss!(
             du, u, p, cache.y, pt, cache.bc, cache.residual, cache.mesh, cache)
@@ -111,48 +113,58 @@ function __construct_nlproblem(cache::MIRKNCache{iip}, y::AbstractVector) where
         @closure (u, p) -> __mirkn_loss(u, p, cache.y, pt, cache.bc, cache.mesh, cache)
     end
 
-    lossₚ = (iip ? __Fix3 : Base.Fix2)(loss, cache.p)
-    sd = alg.jac_alg.diffmode isa AutoSparse ? SymbolicsSparsityDetection() :
-         NoSparsityDetection()
-    ad = alg.jac_alg.diffmode
-    lz = reduce(vcat, cache.y₀)
-    jac_cache = __sparse_jacobian_cache(Val(iip), ad, sd, lossₚ, lz, lz)
-    jac_prototype = init_jacobian(jac_cache)
+    diffmode = if alg.jac_alg.diffmode isa AutoSparse
+        AutoSparse(get_dense_ad(alg.jac_alg.diffmode);
+            # Local because iszero require primal values
+            sparsity_detector = SparseConnectivityTracer.TracerLocalSparsityDetector(), coloring_algorithm = GreedyColoringAlgorithm())
+    else
+        alg.jac_alg.diffmode
+    end
+
+    resid_prototype = similar(y)
+
+    jac_cache = if iip
+        DI.prepare_jacobian(loss, resid_prototype, diffmode, y, Constant(cache.p))
+    else
+        DI.prepare_jacobian(loss, diffmode, y, Constant(cache.p))
+    end
+
+    jac_prototype = if iip
+        DI.jacobian(loss, resid_prototype, jac_cache, diffmode, y, Constant(cache.p))
+    else
+        DI.jacobian(loss, jac_cache, diffmode, y, Constant(cache.p))
+    end
+
     jac = if iip
-        @closure (J, u, p) -> __mirkn_mpoint_jacobian!(J, u, ad, jac_cache, lossₚ, lz)
+        @closure (J, u, p) -> __mirkn_mpoint_jacobian!(
+            J, u, diffmode, jac_cache, loss, resid_prototype, cache.p)
     else
-        @closure (u, p) -> __mirkn_mpoint_jacobian(jac_prototype, u, ad, jac_cache, lossₚ)
+        @closure (u, p) -> __mirkn_mpoint_jacobian(
+            jac_prototype, u, diffmode, jac_cache, loss, cache.p)
     end
-    resid_prototype = zero(lz)
-    _nlf = NonlinearFunction{iip}(
+
+    nlf = NonlinearFunction{iip}(
         loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    nlprob::NonlinearProblem = NonlinearProblem(_nlf, lz, cache.p)
-    return nlprob
+    __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
 end
 
-function __mirkn_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid) where {L}
-    sparse_jacobian!(J, diffmode, diffcache, loss_fn, resid, x)
+function __mirkn_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid, p) where {L}
+    DI.jacobian!(loss_fn, resid, J, diffcache, diffmode, x, Constant(p))
     return J
 end
 
-function __mirkn_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L) where {L}
-    sparse_jacobian!(J, diffmode, diffcache, loss_fn, x)
+function __mirkn_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L, p) where {L}
+    DI.jacobian!(loss_fn, J, diffcache, diffmode, x, Constant(p))
     return J
 end
 
-@inline function __internal_nlsolve_problem(
-        ::SecondOrderBVProblem{uType, tType, iip, nlls}, resid_prototype,
-        u0, args...; kwargs...) where {uType, tType, iip, nlls}
-    return NonlinearProblem(args...; kwargs...)
-end
-
-function __mirkn_mpoint_jacobian!(J, x, diffmode, diffcache, loss, resid)
-    sparse_jacobian!(J, diffmode, diffcache, loss, resid, x)
+function __mirkn_mpoint_jacobian!(J, x, diffmode, diffcache, loss, resid, p)
+    DI.jacobian!(loss, resid, J, diffcache, diffmode, x, Constant(p))
     return nothing
 end
 
-function __mirkn_mpoint_jacobian(J, x, diffmode, diffcache, loss)
-    sparse_jacobian!(J, diffmode, diffcache, loss, x)
+function __mirkn_mpoint_jacobian(J, x, diffmode, diffcache, loss, p)
+    DI.jacobian!(loss, J, diffcache, diffmode, x, Constant(p))
     return J
 end
 
@@ -186,8 +198,9 @@ end
         bc!::BC, residual, mesh, cache::MIRKNCache) where {BC}
     y_ = recursive_unflatten!(y, u)
     resids = [get_tmp(r, u) for r in residual]
+    soly_ = VectorOfArray(y_)
     Φ!(resids[3:end], cache, y_, u, p)
-    eval_bc_residual!(resids, pt, bc!, y_, p, mesh)
+    eval_bc_residual!(resids[1:2], pt, bc!, soly_, p, mesh)
     recursive_flatten!(resid, resids)
     return nothing
 end
@@ -195,7 +208,8 @@ end
 @views function __mirkn_loss(u, p, y, pt::TwoPointSecondOrderBVProblem,
         bc!::BC, mesh, cache::MIRKNCache) where {BC}
     y_ = recursive_unflatten!(y, u)
+    soly_ = VectorOfArray(y_)
     resid_co = Φ(cache, y_, u, p)
-    resid_bc = eval_bc_residual(pt, bc!, y_, p, mesh)
+    resid_bc = eval_bc_residual(pt, bc!, soly_, p, mesh)
     return vcat(resid_bc, mapreduce(vec, vcat, resid_co))
 end
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/utils.jl b/lib/BoundaryValueDiffEqMIRKN/src/utils.jl
index 809d29b1..27abcf71 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/utils.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/utils.jl
@@ -1,96 +1,3 @@
-function general_eval_bc_residual(
-        ::StandardSecondOrderBVProblem, bc::BC, y, p, mesh) where {BC}
-    M = length(y[1])
-    L = length(mesh)
-    res_bc = bc(y[(L + 1):end], y[1:L], p, mesh)
-    return res_bc
-end
-function eval_bc_residual(::StandardSecondOrderBVProblem, bc::BC, y, dy, p, mesh) where {BC}
-    res_bc = bc(dy, y, p, mesh)
-    return res_bc
-end
-function eval_bc_residual(::TwoPointSecondOrderBVProblem, bc::BC, sol, p, mesh) where {BC}
-    M = length(sol[1])
-    L = length(mesh)
-    ua = sol isa AbstractVector ? sol[1] : sol(first(t))[1:M]
-    ub = sol isa AbstractVector ? sol[L] : sol(last(t))[1:M]
-    dua = sol isa AbstractVector ? sol[L + 1] : sol(first(t))[(M + 1):end]
-    dub = sol isa AbstractVector ? sol[end] : sol(last(t))[(M + 1):end]
-    return vcat(bc[1](dua, ua, p), bc[2](dub, ub, p))
-end
-
-function general_eval_bc_residual!(
-        resid, ::StandardSecondOrderBVProblem, bc!::BC, y, p, mesh) where {BC}
-    L = length(mesh)
-    bc!(resid, y[(L + 1):end], y[1:L], p, mesh)
-end
-function eval_bc_residual!(
-        resid, ::StandardSecondOrderBVProblem, bc!::BC, y, dy, p, mesh) where {BC}
-    M = length(resid[1])
-    res_bc = vcat(resid[1], resid[2])
-    bc!(res_bc, dy, y, p, mesh)
-    copyto!(resid[1], res_bc[1:M])
-    copyto!(resid[2], res_bc[(M + 1):end])
-end
-function eval_bc_residual!(
-        resid, ::TwoPointSecondOrderBVProblem, bc::BC, sol, p, mesh) where {BC}
-    M = length(sol[1])
-    L = length(mesh)
-    ua = sol isa AbstractVector ? sol[1] : sol(first(t))[1:M]
-    ub = sol isa AbstractVector ? sol[L] : sol(last(t))[1:M]
-    dua = sol isa AbstractVector ? sol[L + 1] : sol(first(t))[(M + 1):end]
-    dub = sol isa AbstractVector ? sol[end] : sol(last(t))[(M + 1):end]
-    bc[1](resid[1], dua, ua, p)
-    bc[2](resid[2], dub, ub, p)
-end
-
-function __get_bcresid_prototype(::TwoPointSecondOrderBVProblem, prob::BVProblem, u)
-    prototype = if prob.f.bcresid_prototype !== nothing
-        prob.f.bcresid_prototype.x
-    else
-        first(prob.f.bc)(u, prob.p), last(prob.f.bc)(u, prob.p)
-    end
-    return prototype, size.(prototype)
-end
-function __get_bcresid_prototype(::StandardSecondOrderBVProblem, prob::BVProblem, u)
-    prototype = prob.f.bcresid_prototype !== nothing ? prob.f.bcresid_prototype :
-                __zeros_like(u)
-    return prototype, size(prototype)
-end
-
-# Restructure Non-Vector Inputs
-function __vec_f!(ddu, du, u, p, t, f!, u_size)
-    f!(reshape(ddu, u_size), reshape(du, u_size), reshape(u, u_size), p, t)
-    return nothing
-end
-
-__vec_f(du, u, p, t, f, u_size) = vec(f(reshape(du, u_size), reshape(u, u_size), p, t))
-
-function __vec_so_bc!(resid, dsol, sol, p, t, bc!, resid_size, u_size)
-    bc!(reshape(resid, resid_size), __restructure_sol(dsol, u_size),
-        __restructure_sol(sol, u_size), p, t)
-    return nothing
-end
-
-function __vec_so_bc!(resid, dsol, sol, p, bc!, resid_size, u_size)
-    bc!(reshape(resid, resid_size), reshape(dsol, u_size), reshape(sol, u_size), p)
-    return nothing
-end
-
-function __vec_so_bc(dsol, sol, p, t, bc, u_size)
-    vec(bc(__restructure_sol(dsol, u_size), __restructure_sol(sol, u_size), p, t))
-end
-function __vec_so_bc(dsol, sol, p, bc, u_size)
-    vec(bc(reshape(dsol, u_size), reshape(sol, u_size), p))
-end
-
-__get_non_sparse_ad(ad::AbstractADType) = ad
-
-@inline function __initial_guess_on_mesh(
-        prob::SecondOrderBVProblem, u₀::AbstractArray, Nig, p, alias_u0::Bool)
-    return [copy(vec(u₀)) for _ in 1:(2 * (Nig + 1))]
-end
-
 function concrete_jacobian_algorithm(
         jac_alg::BVPJacobianAlgorithm, prob::AbstractBVProblem, alg)
     return concrete_jacobian_algorithm(jac_alg, prob.problem_type, prob, alg)
diff --git a/lib/BoundaryValueDiffEqShooting/Project.toml b/lib/BoundaryValueDiffEqShooting/Project.toml
index ad454829..51182d46 100644
--- a/lib/BoundaryValueDiffEqShooting/Project.toml
+++ b/lib/BoundaryValueDiffEqShooting/Project.toml
@@ -10,12 +10,12 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -25,7 +25,8 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 
 [compat]
 ADTypes = "1.9"
@@ -37,6 +38,7 @@ BoundaryValueDiffEqCore = "1"
 ConcreteStructs = "0.2.3"
 DiffEqBase = "6.158.3"
 DiffEqDevTools = "2.44"
+DifferentiationInterface = "0.6.24"
 FastAlmostBandedMatrices = "0.1.4"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.38"
@@ -45,7 +47,6 @@ InteractiveUtils = "<0.0.1, 1"
 JET = "0.9.12"
 LinearAlgebra = "1.10"
 LinearSolve = "2.36.2"
-Logging = "1.10"
 ODEInterface = "0.5"
 OrdinaryDiffEq = "6.90.1"
 Pkg = "1.10.0"
@@ -59,7 +60,8 @@ Reexport = "1.2"
 SciMLBase = "2.59.1"
 Setfield = "1.1.1"
 SparseArrays = "1.10"
-SparseDiffTools = "2.23"
+SparseConnectivityTracer = "0.6.9"
+SparseMatrixColorings = "0.4.10"
 StaticArrays = "1.9.8"
 Test = "1.10"
 julia = "1.10"
diff --git a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
index e9f8f3b0..1f5697a8 100644
--- a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
@@ -1,47 +1,51 @@
 module BoundaryValueDiffEqShooting
 
-import PrecompileTools: @compile_workload, @setup_workload
-
-using ADTypes, Adapt, ArrayInterface, BoundaryValueDiffEqCore, DiffEqBase, ForwardDiff,
-      LinearAlgebra, Preferences, RecursiveArrayTools, Reexport, SciMLBase, Setfield,
-      SparseDiffTools
-
-using PreallocationTools: PreallocationTools, DiffCache
-
-# Special Matrix Types
-using BandedMatrices, FastAlmostBandedMatrices, SparseArrays
-
-import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
-                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                                __concrete_nonlinearsolve_algorithm, diff!, __any_sparse_ad,
-                                __FastShortcutBVPCompatibleNonlinearPolyalg,
-                                __FastShortcutBVPCompatibleNLLSPolyalg, __cache_trait,
-                                concrete_jacobian_algorithm, eval_bc_residual,
-                                eval_bc_residual!, get_tmp, __maybe_matmul!,
-                                __append_similar!, __extract_problem_details,
-                                __initial_guess, __default_nonsparse_ad,
-                                __maybe_allocate_diffcache, __get_bcresid_prototype,
-                                __similar, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
-                                __materialize_jacobian_algorithm,
-                                recursive_flatten_twopoint!, __internal_nlsolve_problem,
-                                NoDiffCacheNeeded, DiffCacheNeeded, __extract_mesh,
-                                __extract_u0, __has_initial_guess, __initial_guess_length,
-                                __initial_guess_on_mesh, __flatten_initial_guess,
-                                __get_non_sparse_ad, __build_solution, __Fix3,
-                                __sparse_jacobian_cache, __sparsity_detection_alg,
-                                _sparse_like, ColoredMatrix
-
-import ADTypes: AbstractADType
-import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
-import ConcreteStructs: @concrete
-using DiffEqBase: solve
-import FastClosures: @closure
-import ForwardDiff: ForwardDiff, pickchunksize
-import Logging
-import RecursiveArrayTools: ArrayPartition, DiffEqArray
-import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
-
-@reexport using ADTypes, BoundaryValueDiffEqCore, OrdinaryDiffEq, SparseDiffTools, SciMLBase
+using ADTypes
+using ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
+using BandedMatrices: BandedMatrix, Ones
+using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
+                               recursive_flatten, recursive_flatten!, recursive_unflatten!,
+                               __concrete_nonlinearsolve_algorithm, diff!, __any_sparse_ad,
+                               __FastShortcutBVPCompatibleNonlinearPolyalg,
+                               __FastShortcutBVPCompatibleNLLSPolyalg, __cache_trait,
+                               concrete_jacobian_algorithm, eval_bc_residual,
+                               eval_bc_residual!, get_tmp, __maybe_matmul!,
+                               __append_similar!, __extract_problem_details,
+                               __initial_guess, __default_nonsparse_ad,
+                               __maybe_allocate_diffcache, __get_bcresid_prototype,
+                               __similar, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
+                               __materialize_jacobian_algorithm,
+                               recursive_flatten_twopoint!, __internal_nlsolve_problem,
+                               NoDiffCacheNeeded, DiffCacheNeeded, __extract_mesh,
+                               __extract_u0, __has_initial_guess, __initial_guess_length,
+                               __initial_guess_on_mesh, __flatten_initial_guess,
+                               __get_non_sparse_ad, __build_solution, __Fix3, _sparse_like,
+                               get_dense_ad
+
+using ConcreteStructs: @concrete
+using DiffEqBase: DiffEqBase, solve
+using DifferentiationInterface: DifferentiationInterface, Constant, prepare_jacobian
+using FastAlmostBandedMatrices: AlmostBandedMatrix, fillpart, exclusive_bandpart,
+                                finish_part_setindex!
+using FastClosures: @closure
+using ForwardDiff: ForwardDiff, pickchunksize
+using LinearAlgebra
+using PrecompileTools: @compile_workload, @setup_workload
+using Preferences: Preferences
+using Reexport: @reexport
+using RecursiveArrayTools: ArrayPartition, DiffEqArray, VectorOfArray
+using SciMLBase: SciMLBase, AbstractDiffEqInterpolation, StandardBVProblem, __solve,
+                 _unwrap_val
+using Setfield: @set!, @set
+using SparseArrays: sparse
+using SparseConnectivityTracer: SparseConnectivityTracer
+using SparseMatrixColorings: ColoringProblem, GreedyColoringAlgorithm, sparsity_pattern,
+                             ConstantColoringAlgorithm, row_colors, column_colors, coloring,
+                             LargestFirst
+
+const DI = DifferentiationInterface
+
+@reexport using ADTypes, BoundaryValueDiffEqCore, SciMLBase
 
 include("algorithms.jl")
 include("single_shooting.jl")
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index 7c69875e..06f8cc37 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -23,6 +23,7 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa
     else
         __alg
     end
+
     nshoots = alg.nshoots
 
     if prob.problem_type isa TwoPointBVProblem
@@ -92,11 +93,6 @@ function __solve_nlproblem!(
         nodes, cur_nshoot::Int, M::Int, N::Int, resida_len::Int, residb_len::Int,
         solve_internal_odes!::S, bca::B1, bcb::B2, prob, u0, ode_cache_loss_fn,
         ensemblealg, internal_ode_kwargs; kwargs...) where {B1, B2, S}
-    if __any_sparse_ad(alg.jac_alg)
-        J_proto = __generate_sparse_jacobian_prototype(
-            alg, prob.problem_type, bcresid_prototype, u0, N, cur_nshoot)
-    end
-
     resid_prototype = vcat(
         bcresid_prototype[1], similar(u_at_nodes, cur_nshoot * N), bcresid_prototype[2])
 
@@ -104,21 +100,34 @@ function __solve_nlproblem!(
         du, u, p, cur_nshoot, nodes, prob, solve_internal_odes!,
         resida_len, residb_len, N, bca, bcb, ode_cache_loss_fn)
 
-    sd_bvp = alg.jac_alg.diffmode isa AutoSparse ? __sparsity_detection_alg(J_proto) :
-             NoSparsityDetection()
+    diffmode = if alg.jac_alg.diffmode isa AutoSparse
+        colored_result = __generate_sparse_jacobian_prototype(
+            alg, prob.problem_type, bcresid_prototype,
+            u0, N, cur_nshoot, alg.jac_alg.diffmode)
+        partition = ifelse(
+            ADTypes.mode(alg.jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode, :row, :column)
+        constant_matrix_coloring = ifelse(
+            ADTypes.mode(alg.jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode,
+            row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(alg.jac_alg.diffmode),
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
+    else
+        alg.jac_alg.diffmode
+    end
 
     resid_prototype_cached = similar(resid_prototype)
-    jac_cache = sparse_jacobian_cache(
-        alg.jac_alg.diffmode, sd_bvp, nothing, resid_prototype_cached, u_at_nodes)
-    jac_prototype = init_jacobian(jac_cache)
+    jac_cache = DI.prepare_jacobian(nothing, resid_prototype_cached, diffmode, u_at_nodes)
 
     ode_cache_jac_fn = __multiple_shooting_init_jacobian_odecache(
-        ensemblealg, prob, jac_cache, __cache_trait(alg.jac_alg.diffmode),
-        alg.ode_alg, cur_nshoot, u0; internal_ode_kwargs...)
+        ensemblealg, prob, jac_cache, diffmode, alg.ode_alg,
+        cur_nshoot, u0; internal_ode_kwargs...)
 
     loss_fnₚ = @closure (du, u) -> __multiple_shooting_2point_loss!(
         du, u, prob.p, cur_nshoot, nodes, prob, solve_internal_odes!,
         resida_len, residb_len, N, bca, bcb, ode_cache_jac_fn)
+    jac_prototype = DI.jacobian(loss_fnₚ, resid_prototype, jac_cache, diffmode, u_at_nodes)
 
     jac_fn = @closure (J, u, p) -> __multiple_shooting_2point_jacobian!(
         J, u, p, jac_cache, loss_fnₚ, resid_prototype_cached, alg)
@@ -139,10 +148,6 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
         u_at_nodes, nodes, cur_nshoot::Int, M::Int, N::Int, resid_len::Int,
         solve_internal_odes!::S, bc::BC, prob, f::F, u0_size, u0, ode_cache_loss_fn,
         ensemblealg, internal_ode_kwargs; kwargs...) where {BC, F, S}
-    if __any_sparse_ad(alg.jac_alg)
-        J_proto = __generate_sparse_jacobian_prototype(
-            alg, prob.problem_type, bcresid_prototype, u0, N, cur_nshoot)
-    end
     resid_prototype = vcat(bcresid_prototype, similar(u_at_nodes, cur_nshoot * N))
 
     __resid_nodes = resid_prototype[(end - cur_nshoot * N + 1):end]
@@ -154,25 +159,42 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
         N, f, bc, u0_size, prob.tspan, alg.ode_alg, u0, ode_cache_loss_fn)
 
     # ODE Part
-    sd_ode = alg.jac_alg.nonbc_diffmode isa AutoSparse ? __sparsity_detection_alg(J_proto) :
-             NoSparsityDetection()
-    ode_jac_cache = sparse_jacobian_cache(alg.jac_alg.nonbc_diffmode, sd_ode, nothing,
-        similar(u_at_nodes, cur_nshoot * N), u_at_nodes)
+    nonbc_diffmode = if alg.jac_alg.nonbc_diffmode isa AutoSparse
+        colored_result = __generate_sparse_jacobian_prototype(
+            alg, prob.problem_type, bcresid_prototype, u0,
+            N, cur_nshoot, alg.jac_alg.nonbc_diffmode)
+        partition = ifelse(
+            ADTypes.mode(alg.jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode, :row, :column)
+        constant_matrix_coloring = ifelse(
+            ADTypes.mode(alg.jac_alg.nonbc_diffmode) isa ADTypes.ReverseMode,
+            row_colors, column_colors)(colored_result)
+        AutoSparse(get_dense_ad(alg.jac_alg.nonbc_diffmode),
+            sparsity_detector = ADTypes.KnownJacobianSparsityDetector(sparsity_pattern(colored_result)),
+            coloring_algorithm = ConstantColoringAlgorithm{partition}(
+                sparsity_pattern(colored_result), constant_matrix_coloring))
+    else
+        alg.jac_alg.nonbc_diffmode
+    end
+    ode_jac_cache = DI.prepare_jacobian(
+        nothing, similar(u_at_nodes, cur_nshoot * N), nonbc_diffmode, u_at_nodes)
     ode_cache_ode_jac_fn = __multiple_shooting_init_jacobian_odecache(
-        ensemblealg, prob, ode_jac_cache, __cache_trait(alg.jac_alg.nonbc_diffmode),
+        ensemblealg, prob, ode_jac_cache, nonbc_diffmode,
         alg.ode_alg, cur_nshoot, u0; internal_ode_kwargs...)
 
     # BC Part
-    sd_bc = alg.jac_alg.bc_diffmode isa AutoSparse ? SymbolicsSparsityDetection() :
-            NoSparsityDetection()
-    bc_jac_cache = sparse_jacobian_cache(
-        alg.jac_alg.bc_diffmode, sd_bc, nothing, similar(bcresid_prototype), u_at_nodes)
+    bc_diffmode = if alg.jac_alg.bc_diffmode isa AutoSparse
+        AutoSparse(get_dense_ad(alg.jac_alg.bc_diffmode),
+            sparsity_detector = SparseConnectivityTracer.TracerSparsityDetector(),
+            coloring_algorithm = GreedyColoringAlgorithm())
+    else
+        alg.jac_alg.bc_diffmode
+    end
+    bc_jac_cache = DI.prepare_jacobian(
+        nothing, similar(bcresid_prototype), bc_diffmode, u_at_nodes)
     ode_cache_bc_jac_fn = __multiple_shooting_init_jacobian_odecache(
-        ensemblealg, prob, bc_jac_cache, __cache_trait(alg.jac_alg.bc_diffmode),
+        ensemblealg, prob, bc_jac_cache, bc_diffmode,
         alg.ode_alg, cur_nshoot, u0; internal_ode_kwargs...)
 
-    jac_prototype = vcat(init_jacobian(bc_jac_cache), init_jacobian(ode_jac_cache))
-
     # Define the functions now
     ode_fn = @closure (du, u) -> solve_internal_odes!(
         du, u, prob.p, cur_nshoot, nodes, ode_cache_ode_jac_fn)
@@ -180,9 +202,15 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
         du, u, prob.p, cur_nshoot, nodes, prob, solve_internal_odes!, N,
         f, bc, u0_size, prob.tspan, alg.ode_alg, u0, ode_cache_bc_jac_fn)
 
+    jac_prototype_ode = DI.jacobian(ode_fn, similar(u_at_nodes, cur_nshoot * N),
+        ode_jac_cache, nonbc_diffmode, u_at_nodes)
+    jac_prototype_bc = DI.jacobian(
+        bc_fn, similar(bcresid_prototype), bc_jac_cache, bc_diffmode, u_at_nodes)
+    jac_prototype = vcat(jac_prototype_ode, jac_prototype_bc)
+
     jac_fn = @closure (J, u, p) -> __multiple_shooting_mpoint_jacobian!(
-        J, u, p, similar(bcresid_prototype), resid_nodes,
-        ode_jac_cache, bc_jac_cache, ode_fn, bc_fn, alg, N, M)
+        J, u, p, similar(bcresid_prototype), resid_nodes, ode_jac_cache,
+        bc_jac_cache, ode_fn, bc_fn, nonbc_diffmode, bc_diffmode, N, M)
 
     loss_function! = NonlinearFunction{true}(loss_fn; resid_prototype = resid_prototype,
         jac_prototype = jac_prototype, jac = jac_fn)
@@ -209,15 +237,25 @@ function __multiple_shooting_init_odecache(
 end
 
 function __multiple_shooting_init_jacobian_odecache(
-        ensemblealg, prob, jac_cache, ::NoDiffCacheNeeded, alg, nshoots, u; kwargs...)
+        ensemblealg, prob, jac_cache, diffmode, alg, nshoots, u; kwargs...)
+    __multiple_shooting_init_jacobian_odecache(
+        ensemblealg, prob, jac_cache, __cache_trait(diffmode),
+        diffmode, alg, nshoots, u; kwargs...)
+end
+
+function __multiple_shooting_init_jacobian_odecache(
+        ensemblealg, prob, jac_cache, ::NoDiffCacheNeeded,
+        diffmode, alg, nshoots, u; kwargs...)
     return __multiple_shooting_init_odecache(ensemblealg, prob, alg, u, nshoots; kwargs...)
 end
 
 function __multiple_shooting_init_jacobian_odecache(
-        ensemblealg, prob, jac_cache, ::DiffCacheNeeded, alg, nshoots, u; kwargs...)
-    cache = jac_cache.cache
-    if cache isa ForwardDiff.JacobianConfig
-        xduals = reshape(cache.duals[2][1:length(u)], size(u))
+        ensemblealg, prob, jac_cache, ::DiffCacheNeeded,
+        diffmode, alg, nshoots, u; kwargs...)
+    if diffmode isa AutoSparse
+        xduals = reshape(jac_cache.pushforward_prep.xdual_tmp[1:length(u)], size(u))
+    elseif diffmode isa AutoForwardDiff
+        xduals = reshape(jac_cache.config.duals[2][1:length(u)], size(u))
     else
         xduals = reshape(cache.t[1:length(u)], size(u))
     end
@@ -278,19 +316,18 @@ end
 
 function __multiple_shooting_2point_jacobian!(
         J, us, p, jac_cache, loss_fn::F, resid, alg::MultipleShooting) where {F}
-    sparse_jacobian!(J, alg.jac_alg.diffmode, jac_cache, loss_fn, resid, us)
+    DI.jacobian!(loss_fn, resid, J, jac_cache, alg.jac_alg.diffmode, us)
     return nothing
 end
 
 function __multiple_shooting_mpoint_jacobian!(
         J, us, p, resid_bc, resid_nodes, ode_jac_cache, bc_jac_cache, ode_fn::F1,
-        bc_fn::F2, alg::MultipleShooting, N::Int, M::Int) where {F1, F2}
+        bc_fn::F2, nonbc_diffmode, bc_diffmode, N::Int, M::Int) where {F1, F2}
     J_bc = @view(J[1:M, :])
     J_c = @view(J[(M + 1):end, :])
 
-    sparse_jacobian!(
-        J_c, alg.jac_alg.nonbc_diffmode, ode_jac_cache, ode_fn, resid_nodes.du, us)
-    sparse_jacobian!(J_bc, alg.jac_alg.bc_diffmode, bc_jac_cache, bc_fn, resid_bc, us)
+    DI.jacobian!(ode_fn, resid_nodes.du, J_c, ode_jac_cache, nonbc_diffmode, us)
+    DI.jacobian!(bc_fn, resid_bc, J_bc, bc_jac_cache, bc_diffmode, us)
 
     return nothing
 end
diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index 231c8633..55264156 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -34,25 +34,24 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;
             u, p, ode_cache_loss_fn, bc, u0_size, prob.problem_type)
     end
 
-    sd = alg.jac_alg.diffmode isa AutoSparse ? SymbolicsSparsityDetection() :
-         NoSparsityDetection()
     y_ = similar(resid_prototype)
 
-    # Construct the jacobian function
-    # NOTE: We pass in a separate Jacobian Function because that allows us to cache the
-    #       the internal ode solve cache. This cache needs to be distinct from the regular
-    #       residual function cache
+    diffmode = if alg.jac_alg.diffmode isa AutoSparse
+        AutoSparse(get_dense_ad(alg.jac_alg.diffmode),
+            sparsity_detector = SparseConnectivityTracer.TracerSparsityDetector(),
+            coloring_algorithm = GreedyColoringAlgorithm())
+    else
+        alg.jac_alg.diffmode
+    end
+
     jac_cache = if iip
-        sparse_jacobian_cache(alg.jac_alg.diffmode, sd, nothing, y_, vec(u0))
+        DI.prepare_jacobian(nothing, y_, diffmode, vec(u0))
     else
-        sparse_jacobian_cache(alg.jac_alg.diffmode, sd, nothing, vec(u0); fx = y_)
+        DI.prepare_jacobian(nothing, diffmode, vec(u0))
     end
 
     ode_cache_jac_fn = __single_shooting_jacobian_ode_cache(
-        internal_prob, jac_cache, __cache_trait(alg.jac_alg.diffmode),
-        u0, alg.ode_alg; ode_kwargs...)
-
-    jac_prototype = init_jacobian(jac_cache)
+        internal_prob, jac_cache, __cache_trait(diffmode), u0, alg.ode_alg; ode_kwargs...)
 
     loss_fnₚ = if iip
         @closure (du, u) -> __single_shooting_loss!(
@@ -62,12 +61,18 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;
             u, prob.p, ode_cache_jac_fn, bc, u0_size, prob.problem_type)
     end
 
+    jac_prototype = if iip
+        DI.jacobian(loss_fnₚ, y_, jac_cache, diffmode, vec(u0))
+    else
+        DI.jacobian(loss_fnₚ, jac_cache, diffmode, vec(u0))
+    end
+
     jac_fn = if iip
         @closure (J, u, p) -> __single_shooting_jacobian!(
-            J, u, jac_cache, alg.jac_alg.diffmode, loss_fnₚ, y_)
+            J, u, jac_cache, diffmode, loss_fnₚ, y_)
     else
         @closure (u, p) -> __single_shooting_jacobian(
-            jac_prototype, u, jac_cache, alg.jac_alg.diffmode, loss_fnₚ)
+            jac_prototype, u, jac_cache, diffmode, loss_fnₚ)
     end
 
     nlf = NonlinearFunction{iip}(loss_fn; jac_prototype = jac_prototype,
@@ -119,12 +124,12 @@ function __single_shooting_loss(u, p, cache, bc::BC, u0_size, pt) where {BC}
 end
 
 function __single_shooting_jacobian!(J, u, jac_cache, diffmode, loss_fn::L, fu) where {L}
-    sparse_jacobian!(J, diffmode, jac_cache, loss_fn, fu, vec(u))
+    DI.jacobian!(loss_fn, fu, J, jac_cache, diffmode, vec(u))
     return J
 end
 
 function __single_shooting_jacobian(J, u, jac_cache, diffmode, loss_fn::L) where {L}
-    sparse_jacobian!(J, diffmode, jac_cache, loss_fn, vec(u))
+    DI.jacobian!(loss_fn, J, jac_cache, diffmode, vec(u))
     return J
 end
 
@@ -135,7 +140,7 @@ end
 
 function __single_shooting_jacobian_ode_cache(
         prob, jac_cache, ::DiffCacheNeeded, u0, ode_alg; kwargs...)
-    cache = jac_cache.cache
+    cache = jac_cache.config
     if cache isa ForwardDiff.JacobianConfig
         xduals = cache.duals isa Tuple ? cache.duals[2] : cache.duals
     else
diff --git a/lib/BoundaryValueDiffEqShooting/src/sparse_jacobians.jl b/lib/BoundaryValueDiffEqShooting/src/sparse_jacobians.jl
index 29e85277..72471b87 100644
--- a/lib/BoundaryValueDiffEqShooting/src/sparse_jacobians.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/sparse_jacobians.jl
@@ -13,18 +13,22 @@ Returns a 3-Tuple:
     Two-Point Problem) else `nothing`.
 """
 function __generate_sparse_jacobian_prototype(::MultipleShooting, ::StandardBVProblem,
-        bcresid_prototype, u0, N::Int, nshoots::Int)
+        bcresid_prototype, u0, N::Int, nshoots::Int, ad)
     fast_scalar_indexing(u0) ||
         error("Sparse Jacobians are only supported for Fast Scalar Index-able Arrays")
     J₁ = nshoots * N
     J₂ = (nshoots + 1) * N
     J = BandedMatrix(Ones{eltype(u0)}(J₁, J₂), (N - 1, N + 1))
 
-    return ColoredMatrix(sparse(J), matrix_colors(J'), matrix_colors(J))
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
+    result = coloring(sparse(J), problem, algo)
+    return result
 end
 
 function __generate_sparse_jacobian_prototype(::MultipleShooting, ::TwoPointBVProblem,
-        bcresid_prototype, u0, N::Int, nshoots::Int)
+        bcresid_prototype, u0, N::Int, nshoots::Int, ad)
     fast_scalar_indexing(u0) ||
         error("Sparse Jacobians are only supported for Fast Scalar Index-able Arrays")
 
@@ -39,6 +43,10 @@ function __generate_sparse_jacobian_prototype(::MultipleShooting, ::TwoPointBVPr
     J = BandedMatrix(Ones{eltype(u0)}(J₁, J₂), (max(L₁, L₂) + N - 1, N + 1))
 
     # for underdetermined systems we don't have banded qr implemented. use sparse
-    J₁ < J₂ && return ColoredMatrix(sparse(J), matrix_colors(J'), matrix_colors(J))
-    return ColoredMatrix(J, matrix_colors(J'), matrix_colors(J))
+    problem = ColoringProblem(;
+        partition = ifelse((ADTypes.mode(ad) isa ADTypes.ReverseMode), :row, :column))
+    algo = GreedyColoringAlgorithm()
+    # for underdetermined systems we don't have banded qr implemented. use sparse
+    J₁ < J₂ && return coloring(sparse(J), problem, algo)
+    return coloring(J, problem, algo)
 end
diff --git a/lib/BoundaryValueDiffEqShooting/test/basic_problems_tests.jl b/lib/BoundaryValueDiffEqShooting/test/basic_problems_tests.jl
index 2398dedb..1a7074ca 100644
--- a/lib/BoundaryValueDiffEqShooting/test/basic_problems_tests.jl
+++ b/lib/BoundaryValueDiffEqShooting/test/basic_problems_tests.jl
@@ -1,5 +1,5 @@
 @testitem "Basic Shooting" begin
-    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve, JET
+    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve, JET, OrdinaryDiffEq
 
     SOLVERS = [Shooting(Tsit5(), NewtonRaphson(; autodiff = AutoFiniteDiff())),
         Shooting(Tsit5(), NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 2))),
@@ -126,7 +126,7 @@
 end
 
 @testitem "Shooting with Complex Values" begin
-    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve
+    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve, OrdinaryDiffEq
 
     SOLVERS = [
         Shooting(Vern7(), NewtonRaphson(; autodiff = AutoFiniteDiff())), Shooting(Vern7()),
@@ -160,7 +160,7 @@ end
 end
 
 @testitem "Flow in a Channel" begin
-    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve
+    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve, OrdinaryDiffEq
 
     function flow_in_a_channel!(du, u, p, t)
         R, P = p
@@ -209,7 +209,7 @@ end
 end
 #FIXME: MultipleShooting fails for large out-of-place BVP systems
 @testitem "Ray Tracing" begin
-    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve
+    using BoundaryValueDiffEqShooting, LinearAlgebra, LinearSolve, OrdinaryDiffEq
 
     @inline v(x, y, z, p) = 1 / (4 + cos(p[1] * x) + sin(p[2] * y) - cos(p[3] * z))
     @inline ux(x, y, z, p) = -p[1] * sin(p[1] * x)
diff --git a/lib/BoundaryValueDiffEqShooting/test/nlls_tests.jl b/lib/BoundaryValueDiffEqShooting/test/nlls_tests.jl
index 13ba4dcf..d492b19a 100644
--- a/lib/BoundaryValueDiffEqShooting/test/nlls_tests.jl
+++ b/lib/BoundaryValueDiffEqShooting/test/nlls_tests.jl
@@ -1,5 +1,5 @@
 @testitem "Overconstrained BVP" begin
-    using BoundaryValueDiffEqShooting, LinearAlgebra, JET
+    using BoundaryValueDiffEqShooting, OrdinaryDiffEq, LinearAlgebra, JET
 
     SOLVERS = [
         Shooting(Tsit5(), NewtonRaphson(),
diff --git a/lib/BoundaryValueDiffEqShooting/test/orbital_tests.jl b/lib/BoundaryValueDiffEqShooting/test/orbital_tests.jl
index 471ba0de..62e26874 100644
--- a/lib/BoundaryValueDiffEqShooting/test/orbital_tests.jl
+++ b/lib/BoundaryValueDiffEqShooting/test/orbital_tests.jl
@@ -1,5 +1,5 @@
 @testitem "Lambert's Problem" begin
-    using BoundaryValueDiffEqShooting, LinearAlgebra
+    using BoundaryValueDiffEqShooting, OrdinaryDiffEq, LinearAlgebra
 
     y0 = [-4.7763169762853989E+06, -3.8386398704441520E+05, -5.3500183933132319E+06,
         -5528.612564911408, 1216.8442360202787, 4845.114446429901]
diff --git a/src/BoundaryValueDiffEq.jl b/src/BoundaryValueDiffEq.jl
index a1db8819..90f192a3 100644
--- a/src/BoundaryValueDiffEq.jl
+++ b/src/BoundaryValueDiffEq.jl
@@ -1,29 +1,20 @@
 module BoundaryValueDiffEq
 
-import PrecompileTools: @compile_workload, @setup_workload
-
-using ADTypes, Adapt, ArrayInterface, BoundaryValueDiffEqCore, BoundaryValueDiffEqFIRK,
-      BoundaryValueDiffEqMIRK, BoundaryValueDiffEqShooting, DiffEqBase, ForwardDiff,
-      LinearAlgebra, Preferences, NonlinearSolveFirstOrder, RecursiveArrayTools, Reexport,
-      SciMLBase, Setfield, SparseDiffTools
-
-using PreallocationTools: PreallocationTools, DiffCache
-
-# Special Matrix Types
-using BandedMatrices, FastAlmostBandedMatrices, SparseArrays
-
-import ADTypes: AbstractADType
-import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
-import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm
-import ConcreteStructs: @concrete
-import DiffEqBase: solve
-import FastClosures: @closure
-import ForwardDiff: ForwardDiff, pickchunksize
-import Logging
-import RecursiveArrayTools: ArrayPartition, DiffEqArray
-import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
-
-@reexport using ADTypes, DiffEqBase, OrdinaryDiffEq, SparseDiffTools, SciMLBase
+using ADTypes
+using ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
+using BoundaryValueDiffEqAscher
+using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm
+using BoundaryValueDiffEqFIRK
+using BoundaryValueDiffEqMIRK
+using BoundaryValueDiffEqMIRKN
+using BoundaryValueDiffEqShooting
+using DiffEqBase: DiffEqBase, solve
+using FastClosures: @closure
+using ForwardDiff: ForwardDiff, pickchunksize
+using Reexport: @reexport
+using SciMLBase
+
+@reexport using ADTypes, SciMLBase
 
 include("extension_algs.jl")
 
diff --git a/test/misc/manifolds_tests.jl b/test/misc/manifolds_tests.jl
index d47a5340..9dbdabd1 100644
--- a/test/misc/manifolds_tests.jl
+++ b/test/misc/manifolds_tests.jl
@@ -1,5 +1,5 @@
 @testitem "Manifolds.jl Integration" begin
-    using LinearAlgebra
+    using LinearAlgebra, OrdinaryDiffEq
 
     struct EmbeddedTorus
         R::Float64
diff --git a/test/misc/non_vector_input_tests.jl b/test/misc/non_vector_input_tests.jl
index 95413aa9..fbae8294 100644
--- a/test/misc/non_vector_input_tests.jl
+++ b/test/misc/non_vector_input_tests.jl
@@ -1,5 +1,5 @@
 @testitem "Non-Vector Inputs" begin
-    using LinearAlgebra, NonlinearSolveFirstOrder
+    using LinearAlgebra, NonlinearSolveFirstOrder, OrdinaryDiffEq
 
     for order in (2, 3, 4, 5, 6)
         s = Symbol("MIRK$(order)")
diff --git a/test/misc/type_stability_tests.jl b/test/misc/type_stability_tests.jl
index 19a82793..ecabc9e6 100644
--- a/test/misc/type_stability_tests.jl
+++ b/test/misc/type_stability_tests.jl
@@ -1,5 +1,5 @@
 @testitem "Type Stability" begin
-    using LinearAlgebra, BoundaryValueDiffEq
+    using LinearAlgebra, BoundaryValueDiffEq, OrdinaryDiffEq
 
     f(u, p, t) = [p[1] * u[1] - p[2] * u[1] * u[2], p[3] * u[1] * u[2] - p[4] * u[2]]
     function f!(du, u, p, t)
diff --git a/test/wrappers/odeinterface_tests.jl b/test/wrappers/odeinterface_tests.jl
index e3af08c1..f2991c6f 100644
--- a/test/wrappers/odeinterface_tests.jl
+++ b/test/wrappers/odeinterface_tests.jl
@@ -45,7 +45,7 @@ end
 
 # Just test that it runs. BVPSOL only works with linearly separable BCs.
 @testitem "BVPSOL" setup=[ODEInterfaceWrapperTestSetup] begin
-    using ODEInterface, RecursiveArrayTools, NonlinearSolveFirstOrder
+    using ODEInterface, RecursiveArrayTools, NonlinearSolveFirstOrder, OrdinaryDiffEq
 
     tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan,
         p; bcresid_prototype = (zeros(1), zeros(1)))