diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index cb1566aa93..8feaee42f1 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -20,7 +20,6 @@ steps:
     matrix:
       setup:
         version:
-          - "1.9"
           - "1"
         group:
           - "Regression_I"
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 1524481770..ab02143c1e 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -25,7 +25,6 @@ jobs:
           - ODEInterfaceRegression
           - Multithreading
         version:
-          - '1.9'
           - '1'
     steps:
       - uses: actions/checkout@v4
diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml
index 61beadf612..ae25ad90da 100644
--- a/.github/workflows/Downstream.yml
+++ b/.github/workflows/Downstream.yml
@@ -14,7 +14,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        julia-version: [1.9, 1]
+        julia-version: [1]
         os: [ubuntu-latest]
         package:
           - {user: SciML, repo: DelayDiffEq.jl, group: Interface}
diff --git a/Project.toml b/Project.toml
index 5b5eb4085c..029cfde0cf 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEq"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <mayingbo5@gmail.com>"]
-version = "6.63.0"
+version = "6.69.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -23,7 +23,6 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
-LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
@@ -63,18 +62,17 @@ LineSearches = "7"
 LinearAlgebra = "1.9"
 LinearSolve = "2.1.10"
 Logging = "1.9"
-LoopVectorization = "0.12"
 MacroTools = "0.5"
 MuladdMacro = "0.2.1"
 NLsolve = "4"
-NonlinearSolve = "3"
+NonlinearSolve = "3.3"
 Polyester = "0.7"
-PreallocationTools = "0.4"
+PreallocationTools = "0.4.15"
 PrecompileTools = "1"
 Preferences = "1.3"
 RecursiveArrayTools = "2.36, 3"
 Reexport = "1.0"
-SciMLBase = "2"
+SciMLBase = "2.17"
 SciMLOperators = "0.3"
 SimpleNonlinearSolve = "1"
 SimpleUnPack = "1"
@@ -83,7 +81,7 @@ SparseDiffTools = "2.3"
 StaticArrayInterface = "1.2"
 StaticArrays = "1.0"
 TruncatedStacktraces = "1.2"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
diff --git a/README.md b/README.md
index b4e9bf8ad4..7b07bdbe0c 100644
--- a/README.md
+++ b/README.md
@@ -15,6 +15,7 @@ ordinary differential equation solvers and utilities. While completely independe
 and usable on its own, users interested in using this
 functionality should check out [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl).
 
+
 ## Installation
 
 Assuming that you already have Julia correctly installed, it suffices to import
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 3d23a1f4bc..65bdf16001 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -18,8 +18,6 @@ using MuladdMacro, SparseArrays, FastClosures
 
 using LinearAlgebra
 
-using LoopVectorization
-
 import StaticArrayInterface
 
 import InteractiveUtils
diff --git a/src/dense/generic_dense.jl b/src/dense/generic_dense.jl
index 22fc4860ba..3683c25bed 100644
--- a/src/dense/generic_dense.jl
+++ b/src/dense/generic_dense.jl
@@ -360,6 +360,7 @@ function ode_interpolation(tvals, id::I, idxs, deriv::D, p,
             i₋ = max(1, _searchsortedlast(ts, t, i₋, tdir > 0))
             i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
         end
+        id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
         i₋₊ref[] = (i₋, i₊)
         dt = ts[i₊] - ts[i₋]
         Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t - ts[i₋]) / dt
@@ -405,6 +406,7 @@ function ode_interpolation!(vals, tvals, id::I, idxs, deriv::D, p,
             i₋ = max(1, _searchsortedlast(ts, t, i₋, tdir > 0))
             i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
         end
+        id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
 
         dt = ts[i₊] - ts[i₋]
         Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t - ts[i₋]) / dt
@@ -477,6 +479,7 @@ function ode_interpolation(tval::Number, id::I, idxs, deriv::D, p,
         i₋ = max(1, _searchsortedlast(ts, tval, 1, tdir > 0))
         i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
     end
+    id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
 
     @inbounds begin
         dt = ts[i₊] - ts[i₋]
@@ -525,6 +528,7 @@ function ode_interpolation!(out, tval::Number, id::I, idxs, deriv::D, p,
         i₋ = max(1, _searchsortedlast(ts, tval, 1, tdir > 0))
         i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
     end
+    id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE)
 
     @inbounds begin
         dt = ts[i₊] - ts[i₋]
@@ -687,8 +691,13 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
-    @inbounds (1 - Θ) * y₀ + Θ * y₁ +
-              differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]))
+    if all(differential_vars)
+        @inbounds (1 - Θ) * y₀ + Θ * y₁ +
+                (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]))
+    else
+        @inbounds (1 - Θ) * y₀ + Θ * y₁ +
+                differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2]))
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,
@@ -755,10 +764,17 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
-    @inbounds (.!differential_vars).*(y₁ - y₀)/dt + differential_vars .*(
-              k[1] +
-              Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
-               Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt)
+    if all(differential_vars)
+        @inbounds  (
+            k[1] +
+            Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+            Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt)
+    else
+        @inbounds (.!differential_vars).*(y₁ - y₀)/dt + differential_vars .*(
+                k[1] +
+                Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+                Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt)
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,
@@ -826,8 +842,13 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
-    @inbounds differential_vars .* (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
-               Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt)
+    if all(differential_vars)
+        @inbounds (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+        Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt)
+    else
+        @inbounds differential_vars .* (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ +
+                Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt)
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,
@@ -887,7 +908,11 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing,
     T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite
     #@.. broadcast=false (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
-    @inbounds differential_vars .* (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt)
+    if all(differential_vars)
+        @inbounds (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt)
+    else
+        @inbounds differential_vars .* (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt)
+    end
 end
 
 @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing,
diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl
index 060423ab81..164e178b24 100644
--- a/src/dense/interpolants.jl
+++ b/src/dense/interpolants.jl
@@ -30,16 +30,28 @@ end
 
 @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
     cache::Union{FunctionMapConstantCache, FunctionMapCache},
-    idxs, T::Type{Val{0}}, differential_vars::Nothing)
+    idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing)
     y₀
 end
 
-@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
     cache::Union{FunctionMapConstantCache, FunctionMapCache},
     idxs, T::Type{Val{0}}, differential_vars::Nothing)
+    y₀[idxs]
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{FunctionMapConstantCache, FunctionMapCache},
+    idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing)
     recursivecopy!(out, y₀)
 end
 
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{FunctionMapConstantCache, FunctionMapCache},
+    idxs, T::Type{Val{0}}, differential_vars::Nothing)
+    @views out[idxs] .= y₀[idxs]
+end
+
 """
 Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 192
 """
diff --git a/src/dense/rosenbrock_interpolants.jl b/src/dense/rosenbrock_interpolants.jl
index 2e73d536bc..11156c0334 100644
--- a/src/dense/rosenbrock_interpolants.jl
+++ b/src/dense/rosenbrock_interpolants.jl
@@ -314,4 +314,70 @@ end
                                     y₀[idxs] + y₁[idxs]) / dt
     out
 end
-#-
+
+# Second Derivative
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache,
+    idxs::Nothing, T::Type{Val{2}}, differential_vars)
+    @inbounds (-2 * k[1] + 2 * k[2] + Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3])) / dt^2
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing,
+    T::Type{Val{2}}, differential_vars)
+    @inbounds @.. broadcast=false (-2 * k[1] + 2 * k[2] +
+                                    Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3]))/dt^2
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{2}}, differential_vars)
+    @.. broadcast=false (-2 * k[1][idxs] + 2 * k[2][idxs] +
+                          Θ * (-6 * k[2][idxs] + 6 * k[3][idxs] - 12 * Θ * k[3][idxs]))/dt^2
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs::Nothing, T::Type{Val{2}}, differential_vars)
+    @.. broadcast=false out= (-2 * k[1] + 2 * k[2] +
+                              Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3])) / dt^2
+    out
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{2}}, differential_vars)
+    @views @.. broadcast=false out= (-2 * k[1][idxs] + 2 * k[2][idxs] +
+                                     Θ *
+                                     (-6 * k[2][idxs] + 6 * k[3][idxs] - 12 * Θ * k[3][idxs])) / dt^2
+    out
+end
+
+# Third Derivative
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache,
+    idxs::Nothing, T::Type{Val{3}}, differential_vars)
+    @inbounds (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3]) / dt^3
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing,
+    T::Type{Val{3}}, differential_vars)
+    @inbounds @.. broadcast=false (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3])/dt^3
+end
+
+@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{3}}, differential_vars)
+    @.. broadcast=false (-6 * k[2][idxs] + 6 * k[3][idxs] - 24 * Θ * k[3][idxs])/dt^3
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs::Nothing, T::Type{Val{3}}, differential_vars)
+    @.. broadcast=false out= (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3]) / dt^3
+    out
+end
+
+@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
+    cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache},
+    idxs, T::Type{Val{3}}, differential_vars)
+    @views @.. broadcast=false out= (-6 * k[2][idxs] + 6 * k[3][idxs] - 24 * Θ * k[3][idxs]) / dt^3
+    out
+end
diff --git a/src/derivative_wrappers.jl b/src/derivative_wrappers.jl
index 75f3484407..0123e2e578 100644
--- a/src/derivative_wrappers.jl
+++ b/src/derivative_wrappers.jl
@@ -6,7 +6,7 @@ const FIRST_AUTODIFF_TGRAD_MESSAGE = """
                                1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
                                   Rosenbrock23(autodiff=false)). More details can be found at
                                   https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
-                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic 
+                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic
                                   differentiation (using tools like PreallocationTools.jl). More details
                                   can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
                                3. Defining analytical Jacobians and time gradients. More details can be
@@ -48,7 +48,7 @@ const FIRST_AUTODIFF_JAC_MESSAGE = """
                                1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
                                   Rosenbrock23(autodiff=false)). More details can befound at
                                   https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
-                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic 
+                               2. Improving the compatibility of `f` with ForwardDiff.jl automatic
                                   differentiation (using tools like PreallocationTools.jl). More details
                                   can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
                                3. Defining analytical Jacobians. More details can be
diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl
index bb71431648..9d7537168b 100644
--- a/src/initialize_dae.jl
+++ b/src/initialize_dae.jl
@@ -141,7 +141,20 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation
         failed = nlsolvefail(nlsolver)
         @.. broadcast=false integrator.u=integrator.uprev + z
     else
-        isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+
+        # _u0 should be non-dual since NonlinearSolve does not differentiate the solver
+        # These non-dual values are thus used to make the caches
+        #_du = DiffEqBase.value.(du)
+        _u0 = DiffEqBase.value.(u0)
+
+        # If not doing auto-diff of the solver, save an allocation
+        if typeof(u0) === typeof(_u0)
+            tmp = get_tmp_cache(integrator)[1]
+        else
+            tmp = copy(_u0)
+        end
+
+        isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u0) !== typeof(_u0)
         if isAD
             chunk = ForwardDiff.pickchunksize(length(tmp))
             _tmp = PreallocationTools.dualcache(tmp, chunk)
@@ -150,10 +163,16 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation
         end
 
         nlequation! = @closure (out, u, p) -> begin
+            TP = DiffEqBase.anyeltypedual(p)
+            if TP <: Dual
+                T = Base.promote_type(eltype(u), TP)
+            else
+                T = eltype(u)
+            end
             update_coefficients!(M, u, p, t)
             # f(u,p,t) + M * (u0 - u)/dt
-            tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp
-            @. tmp = (u0 - u) / dt
+            tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp
+            @. tmp = (_u0 - u) / dt
             mul!(_vec(out), M, _vec(tmp))
             f(tmp, u, p, t)
             out .+= tmp
@@ -277,7 +296,6 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     u0 = integrator.u
 
     dtmax = integrator.opts.dtmax
-    tmp = get_tmp_cache(integrator)[1]
     resid = get_tmp_cache(integrator)[2]
 
     dt = t != 0 ? min(t / 1000, dtmax / 10) : dtmax / 10 # Haven't implemented norm reduction
@@ -285,7 +303,19 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     f(resid, integrator.du, u0, p, t)
     integrator.opts.internalnorm(resid, t) <= integrator.opts.abstol && return
 
-    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+    # _du and _u should be non-dual since NonlinearSolve does not differentiate the solver
+    # These non-dual values are thus used to make the caches
+    #_du = DiffEqBase.value.(du)
+    _u0 = DiffEqBase.value.(u0)
+
+    # If not doing auto-diff of the solver, save an allocation
+    if typeof(u0) === typeof(_u0)
+        tmp = get_tmp_cache(integrator)[1]
+    else
+        tmp = copy(_u0)
+    end
+
+    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u0) !== typeof(_u0)
     if isAD
         chunk = ForwardDiff.pickchunksize(length(tmp))
         _tmp = PreallocationTools.dualcache(tmp, chunk)
@@ -294,9 +324,15 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     end
 
     nlequation! = @closure (out, u, p) -> begin
-        tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp
+        TP = DiffEqBase.anyeltypedual(p)
+        if TP <: Dual
+            T = Base.promote_type(eltype(u), TP)
+        else
+            T = eltype(u)
+        end
+        tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp
         #M * (u-u0)/dt - f(u,p,t)
-        @. tmp = (u - u0) / dt
+        @. tmp = (u - _u0) / dt
         f(out, tmp, u, p, t)
         nothing
     end
@@ -414,9 +450,24 @@ function _initialize_dae!(integrator, prob::ODEProblem,
     integrator.opts.internalnorm(tmp, t) <= alg.abstol && return
     alg_u = @view u[algebraic_vars]
 
-    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+    # These non-dual values are thus used to make the caches
+    _u = DiffEqBase.value.(u)
+
+    # If auto-diff of the solver, should be non-dual since NonlinearSolve does not differentiate the solver
+    if typeof(u) !== typeof(_u)
+        tmp = DiffEqBase.value.(tmp)
+    end
+
+    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u) !== typeof(_u)
     if isAD
-        chunk = ForwardDiff.pickchunksize(count(algebraic_vars))
+        csize = count(algebraic_vars)
+        if !(p isa SciMLBase.NullParameters) && typeof(_u) !== typeof(u)
+            try
+                csize = max(csize,length(p))
+            catch
+            end
+        end
+        chunk = ForwardDiff.pickchunksize(csize)
         _tmp = PreallocationTools.dualcache(tmp, chunk)
         _du_tmp = PreallocationTools.dualcache(similar(tmp), chunk)
     else
@@ -424,9 +475,15 @@ function _initialize_dae!(integrator, prob::ODEProblem,
     end
 
     nlequation! = @closure (out, x, p) -> begin
-        uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp
-        du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp
-        copyto!(uu, integrator.u)
+        TP = DiffEqBase.anyeltypedual(p)
+        if TP <: Dual
+            T = Base.promote_type(eltype(x), TP)
+        else
+            T = eltype(x)
+        end
+        uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp
+        du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp
+        copyto!(uu, _u)
         alg_uu = @view uu[algebraic_vars]
         alg_uu .= x
         f(du_tmp, uu, p, t)
@@ -531,31 +588,51 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     u = integrator.u
     du = integrator.du
 
-    tmp = get_tmp_cache(integrator)[1]
-    du_tmp = get_tmp_cache(integrator)[2]
-    f(tmp, du, u, p, t)
+    # _du and _u should be non-dual since NonlinearSolve does not differentiate the solver
+    # These non-dual values are thus used to make the caches
+    _du = DiffEqBase.value.(du)
+    _u = DiffEqBase.value.(u)
+
+    # If not doing auto-diff of the solver, save an allocation
+    if typeof(u) === typeof(_u)
+        tmp = get_tmp_cache(integrator)[1]
+        du_tmp = get_tmp_cache(integrator)[2]
+    else
+        tmp = copy(_u)
+        du_tmp = copy(_du)
+    end
+
+    # Can be the same as tmp
+    normtmp = get_tmp_cache(integrator)[1]
+    f(normtmp, du, u, p, t)
 
-    if integrator.opts.internalnorm(tmp, t) <= alg.abstol
+    if integrator.opts.internalnorm(normtmp, t) <= alg.abstol
         return
     elseif differential_vars === nothing
         error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.")
     end
 
-    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
-    if isAD
-        chunk = ForwardDiff.pickchunksize(length(tmp))
-        _tmp = PreallocationTools.dualcache(tmp, chunk)
-        _du_tmp = PreallocationTools.dualcache(du_tmp, chunk)
-    else
-        _tmp, _du_tmp = tmp, similar(tmp)
-    end
+    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u) !== typeof(_u)
+   if isAD
+       chunk = ForwardDiff.pickchunksize(length(tmp))
+       _tmp = PreallocationTools.dualcache(tmp, chunk)
+       _du_tmp = PreallocationTools.dualcache(du_tmp, chunk)
+   else
+        _tmp, _du_tmp = tmp, _du_tmp
+   end
 
     nlequation! = @closure (out, x, p) -> begin
-        du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp
-        uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp
+        TP = DiffEqBase.anyeltypedual(p)
+        if TP <: Dual
+            T = Base.promote_type(eltype(x), TP)
+        else
+            T = eltype(x)
+        end
+        du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp
+        uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp
 
-        @. du_tmp = ifelse(differential_vars, x, du)
-        @. uu = ifelse(differential_vars, u, x)
+        @. du_tmp = ifelse(differential_vars, x, _du)
+        @. uu = ifelse(differential_vars, _u, x)
 
         f(out, du_tmp, uu, p, t)
     end
diff --git a/src/integrators/integrator_utils.jl b/src/integrators/integrator_utils.jl
index 4bc032d1a6..8b246122b7 100644
--- a/src/integrators/integrator_utils.jl
+++ b/src/integrators/integrator_utils.jl
@@ -79,6 +79,10 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
                     integrator.cache.current)
             end
         else # ==t, just save
+            if curt == integrator.sol.prob.tspan[2] && !integrator.opts.save_end
+                integrator.saveiter -= 1
+                continue
+            end
             savedexactly = true
             copyat_or_push!(integrator.sol.t, integrator.saveiter, integrator.t)
             if integrator.opts.save_idxs === nothing
@@ -107,7 +111,9 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
         end
     end
     if force_save || (integrator.opts.save_everystep &&
-        (isempty(integrator.sol.t) || (integrator.t !== integrator.sol.t[end])))
+        (isempty(integrator.sol.t) || (integrator.t !== integrator.sol.t[end]) &&
+        (integrator.opts.save_end || integrator.t !== integrator.sol.prob.tspan[2])
+        ))
         integrator.saveiter += 1
         saved, savedexactly = true, true
         if integrator.opts.save_idxs === nothing
diff --git a/src/integrators/type.jl b/src/integrators/type.jl
index 3229ff80df..1300ac956a 100644
--- a/src/integrators/type.jl
+++ b/src/integrators/type.jl
@@ -186,12 +186,3 @@ end
 # When this is changed, DelayDiffEq.jl must be changed as well!
 
 TruncatedStacktraces.@truncate_stacktrace ODEIntegrator 2 1 3 4
-
-function Base.getproperty(integ::ODEIntegrator, s::Symbol)
-    if s === :destats
-        @warn "destats has been deprecated for stats"
-        getfield(integ, :stats)
-    else
-        getfield(integ, s)
-    end
-end
diff --git a/src/interp_func.jl b/src/interp_func.jl
index f1a3ed4a4e..fe0279966c 100644
--- a/src/interp_func.jl
+++ b/src/interp_func.jl
@@ -1,27 +1,25 @@
 abstract type OrdinaryDiffEqInterpolation{cacheType} <:
               DiffEqBase.AbstractDiffEqInterpolation end
 
-struct InterpolationData{F, uType, tType, kType, cacheType, DV} <:
+struct InterpolationData{F, uType, tType, kType, algType <: Union{Nothing, Vector{Int}}, cacheType, DV} <:
        OrdinaryDiffEqInterpolation{cacheType}
     f::F
     timeseries::uType
     ts::tType
     ks::kType
+    alg_choice::algType
     dense::Bool
     cache::cacheType
     differential_vars::DV
+    sensitivitymode::Bool
 end
 
-struct CompositeInterpolationData{F, uType, tType, kType, cacheType, DV} <:
-       OrdinaryDiffEqInterpolation{cacheType}
-    f::F
-    timeseries::uType
-    ts::tType
-    ks::kType
-    alg_choice::Vector{Int}
-    dense::Bool
-    cache::cacheType
-    differential_vars::DV
+@static if isdefined(SciMLBase, :enable_interpolation_sensitivitymode)
+    function SciMLBase.enable_interpolation_sensitivitymode(interp::InterpolationData)
+        InterpolationData(interp.f,interp.timeseries,interp.ts,interp.ks,
+                        interp.alg_choice, interp.dense, interp.cache,
+                        interp.differential_vars, true)
+    end
 end
 
 function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{
@@ -169,31 +167,17 @@ end
 function (interp::InterpolationData)(tvals, idxs, deriv, p, continuity::Symbol = :left)
     ode_interpolation(tvals, interp, idxs, deriv, p, continuity)
 end
-function (interp::CompositeInterpolationData)(tvals, idxs, deriv, p,
-    continuity::Symbol = :left)
-    ode_interpolation(tvals, interp, idxs, deriv, p, continuity)
-end
 function (interp::InterpolationData)(val, tvals, idxs, deriv, p, continuity::Symbol = :left)
     ode_interpolation!(val, tvals, interp, idxs, deriv, p, continuity)
 end
-function (interp::CompositeInterpolationData)(val, tvals, idxs, deriv, p,
-    continuity::Symbol = :left)
-    ode_interpolation!(val, tvals, interp, idxs, deriv, p, continuity)
-end
 
 function InterpolationData(id::InterpolationData, f)
     InterpolationData(f, id.timeseries,
-        id.ts,
-        id.ks,
-        id.dense,
-        id.cache)
-end
-
-function CompositeInterpolationData(id::CompositeInterpolationData, f)
-    CompositeInterpolationData(f, id.timeseries,
         id.ts,
         id.ks,
         id.alg_choice,
         id.dense,
-        id.cache)
+        id.cache,
+        id.differential_vars,
+        id.sensitivitymode)
 end
diff --git a/src/nlsolve/nlsolve.jl b/src/nlsolve/nlsolve.jl
index b522239ae3..140fbcec7d 100644
--- a/src/nlsolve/nlsolve.jl
+++ b/src/nlsolve/nlsolve.jl
@@ -89,7 +89,7 @@ function nlsolve!(nlsolver::AbstractNLSolver, integrator::DiffEqBase.DEIntegrato
         apply_step!(nlsolver, integrator)
 
         # check for convergence
-        iter > 1 && (η = θ / (1 - θ))
+        iter > 1 && (η = DiffEqBase.value(θ / (1 - θ)))
         if (iter == 1 && ndz < 1e-5) || (iter > 1 && (η >= zero(η) && η * ndz < κ))
             nlsolver.status = Convergence
             nlsolver.nfails = 0
diff --git a/src/perform_step/symplectic_perform_step.jl b/src/perform_step/symplectic_perform_step.jl
index 44a16bb110..3e4f6f849f 100644
--- a/src/perform_step/symplectic_perform_step.jl
+++ b/src/perform_step/symplectic_perform_step.jl
@@ -568,7 +568,7 @@ end
     kdu = f.f1(du, u, p, tnew)
     du = du + dt * a5 * kdu
 
-    tnew = tnew + t + a5 * dt
+    tnew = tnew + a5 * dt
     ku = f.f2(du, u, p, tnew)
     u = u + dt * b6 * ku
 
@@ -625,7 +625,7 @@ end
     f.f1(kdu, du, u, p, tnew)
     @.. broadcast=false du=du + dt * a5 * kdu
 
-    tnew = tnew + t + a5 * dt
+    tnew = tnew + a5 * dt
     f.f2(ku, du, u, p, tnew)
     @.. broadcast=false u=u + dt * b6 * ku
 
diff --git a/src/solve.jl b/src/solve.jl
index acb6ffee54..c999c84f66 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -31,7 +31,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
                 (dense) || !isempty(saveat), # and no dense output
     dt = alg isa FunctionMap && isempty(tstops) ?
          eltype(prob.tspan)(1) : eltype(prob.tspan)(0),
-    dtmin = nothing,
+    dtmin = eltype(prob.tspan)(0),
     dtmax = eltype(prob.tspan)((prob.tspan[end] - prob.tspan[1])),
     force_dtmin = false,
     adaptive = anyadaptive(alg),
@@ -195,10 +195,10 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
         reltol_internal = false
     elseif reltol === nothing
         if uBottomEltypeNoUnits == uBottomEltype
-            reltol_internal = real(convert(uBottomEltype,
-                oneunit(uBottomEltype) * 1 // 10^3))
+            reltol_internal = ForwardDiff.value(real(convert(uBottomEltype,
+                oneunit(uBottomEltype) * 1 // 10^3)))
         else
-            reltol_internal = real.(oneunit.(u) .* 1 // 10^3)
+            reltol_internal = ForwardDiff.value.(real.(oneunit.(u) .* 1 // 10^3))
         end
     else
         reltol_internal = real.(reltol)
@@ -272,7 +272,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
 
     ts = ts_init === () ? tType[] : convert(Vector{tType}, ts_init)
     ks = ks_init === () ? ksEltype[] : convert(Vector{ksEltype}, ks_init)
-    alg_choice = _alg isa CompositeAlgorithm ? Int[] : ()
+    alg_choice = _alg isa CompositeAlgorithm ? Int[] : nothing
 
     if (!adaptive || !isadaptive(_alg)) && save_everystep && tspan[2] - tspan[1] != Inf
         if dt == 0
@@ -281,7 +281,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
             # For fixed dt, the only time dtmin makes sense is if it's smaller than eps().
             # Therefore user specified dtmin doesn't matter, but we need to ensure dt>=eps()
             # to prevent infinite loops.
-            abs(dt) < DiffEqBase.prob2dtmin(prob) &&
+            abs(dt) < dtmin &&
                 throw(ArgumentError("Supplied dt is smaller than dtmin"))
             steps = ceil(Int, internalnorm((tspan[2] - tspan[1]) / dt, tspan[1]))
         end
@@ -293,9 +293,16 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
         sizehint!(ts, 50)
         sizehint!(ks, 50)
     elseif !isempty(saveat_internal)
-        sizehint!(timeseries, length(saveat_internal) + 1)
-        sizehint!(ts, length(saveat_internal) + 1)
-        sizehint!(ks, length(saveat_internal) + 1)
+        savelength = length(saveat_internal) + 1
+        if save_start == false
+            savelength -= 1
+        end
+        if save_end == false && prob.tspan[2] in saveat_internal.valtree
+            savelength -= 1
+        end
+        sizehint!(timeseries, savelength)
+        sizehint!(ts, savelength)
+        sizehint!(ks, savelength)
     else
         sizehint!(timeseries, 2)
         sizehint!(ts, 2)
@@ -351,8 +358,6 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
         controller = default_controller(_alg, cache, qoldinit, beta1, beta2)
     end
 
-    dtmin === nothing && (dtmin = DiffEqBase.prob2dtmin(prob))
-
     save_end_user = save_end
     save_end = save_end === nothing ?
                save_everystep || isempty(saveat) || saveat isa Number ||
@@ -405,19 +410,12 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
 
     stats = SciMLBase.DEStats(0)
     differential_vars = prob isa DAEProblem ? prob.differential_vars : get_differential_vars(f, u)
-      
-    if _alg isa OrdinaryDiffEqCompositeAlgorithm
-        id = CompositeInterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars)
-        sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries,
-            dense = dense, k = ks, interp = id,
-            alg_choice = alg_choice,
-            calculate_error = false, stats = stats)
-    else
-        id = InterpolationData(f, timeseries, ts, ks, dense, cache, differential_vars)
-        sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries,
-            dense = dense, k = ks, interp = id,
-            calculate_error = false, stats = stats)
-    end
+
+    id = InterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars, false)
+    sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries,
+        dense = dense, k = ks, interp = id,
+        alg_choice = alg_choice,
+        calculate_error = false, stats = stats)
 
     if recompile_flag == true
         FType = typeof(f)
diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl
index 3d496ee4c8..c0442ca591 100644
--- a/test/algconvergence/symplectic_tests.jl
+++ b/test/algconvergence/symplectic_tests.jl
@@ -1,69 +1,99 @@
-
-using Test, LinearAlgebra
-using OrdinaryDiffEq, DiffEqBase
-
-# algorithm, dq(p) != p, convergence order
-const ALGOS = ((SymplecticEuler, true, 1),
-    (VelocityVerlet, false, 2),
-    (VerletLeapfrog, true, 2),
-    (PseudoVerletLeapfrog, true, 2),
-    (McAte2, true, 2),
-    (Ruth3, true, 3),
-    (McAte3, true, 3),
-    (CandyRoz4, true, 4),
-    (McAte4, true, 4),
-    (CalvoSanz4, true, 4),
-    (McAte42, true, 1), # known to be broken
-    (McAte5, true, 5),
-    (Yoshida6, true, 6),
-    (KahanLi6, true, 6),
-    (McAte8, true, 8),
-    (KahanLi8, true, 8),
-    (SofSpa10, true, 10))
-
-function dp(p, q, pa, t)
-    0q .+ pa[2]
-end
-
-function dq(p, q, pa, t)
-    p .* pa[1]
-end
-
-dp(res, p, q, pa, t) = (res .= dp(p, q, pa, t))
-dq(res, p, q, pa, t) = (res .= dq(p, q, pa, t))
-
-dynode(iip, dp, dq) = DynamicalODEFunction{iip}(dp, dq)
-
-# [0:1] used in dp, dq; [3:4] start values for p0, q0
-const PARAMS = ((1.0, 0.1, 1.0, 0.0), (0.1, 1.0, 1.0, -1.0))
-const IIPS = (true, false)
-const TSPAN = (0.0, 1.0)
-
-solution(t, w) = (w[2] * t + w[3], (w[2] / 2 * t + w[3]) * w[1] * t + w[4])
-apa(iip::Bool, x) = iip ? vcat.(x) : x
-errorbound(dt, d, x) = 100 * abs(dt)^d + 1000 * eps(norm(x))
-function printerrors(text, calc, solution, pa, t1)
-    print(text, ": ")
-    print(norm(calc[1] - solution(t1, pa)[1]), " ")
-    print(norm(calc[2] - solution(t1, pa)[2]))
-    println()
-end
-
-@testset "symplectic $alg-$iip-$pa" for (alg, x, d) in ALGOS, iip in IIPS, pa in PARAMS
-    dt = 0.01
-    tspan = TSPAN
-    t0, t1 = tspan
-    dynfun = dynode(iip, dp, dq)
-    p0, q0 = apa(iip, solution(t0, pa))
-    prob = DynamicalODEProblem(dynfun, p0, q0, tspan, pa)
-
-    if x || pa[1] == 1
-        sol = solve(prob, alg(); dt = dt)
-        calc = sol(t1)
-        # printerrors("$alg-$iip-$pa", calc, solution, pa, t1)
-        @test calc[1]≈solution(t1, pa)[1] rtol=errorbound(dt, d, calc[1])
-        @test calc[2]≈solution(t1, pa)[2] rtol=errorbound(dt, d, calc[2])
-    else
-        @test_throws ArgumentError solve(prob, alg(); dt = dt)
-    end
-end
+
+using Test, LinearAlgebra
+using OrdinaryDiffEq, DiffEqBase
+
+# algorithm, dq(p) != p, convergence order
+const ALGOS = ((SymplecticEuler, true, 1),
+    (VelocityVerlet, false, 2),
+    (VerletLeapfrog, true, 2),
+    (PseudoVerletLeapfrog, true, 2),
+    (McAte2, true, 2),
+    (Ruth3, true, 3),
+    (McAte3, true, 3),
+    (CandyRoz4, true, 4),
+    (McAte4, true, 4),
+    (CalvoSanz4, true, 4),
+    (McAte42, true, 1), # known to be broken
+    (McAte5, true, 5),
+    (Yoshida6, true, 6),
+    (KahanLi6, true, 6),
+    (McAte8, true, 8),
+    (KahanLi8, true, 8),
+    (SofSpa10, true, 10))
+
+function dp(p, q, pa, t)
+    0q .+ pa[2]
+end
+
+function dq(p, q, pa, t)
+    p .* pa[1]
+end
+
+dp(res, p, q, pa, t) = (res .= dp(p, q, pa, t))
+dq(res, p, q, pa, t) = (res .= dq(p, q, pa, t))
+
+dynode(iip, dp, dq) = DynamicalODEFunction{iip}(dp, dq)
+
+# [0:1] used in dp, dq; [3:4] start values for p0, q0
+const PARAMS = ((1.0, 0.1, 1.0, 0.0), (0.1, 1.0, 1.0, -1.0))
+const IIPS = (true, false)
+const TSPAN = (0.0, 1.0)
+
+solution(t, w) = (w[2] * t + w[3], (w[2] / 2 * t + w[3]) * w[1] * t + w[4])
+apa(iip::Bool, x) = iip ? vcat.(x) : x
+errorbound(dt, d, x) = 100 * abs(dt)^d + 1000 * eps(norm(x))
+function printerrors(text, calc, solution, pa, t1)
+    print(text, ": ")
+    print(norm(calc[1] - solution(t1, pa)[1]), " ")
+    print(norm(calc[2] - solution(t1, pa)[2]))
+    println()
+end
+
+@testset "symplectic $alg-$iip-$pa" for (alg, x, d) in ALGOS, iip in IIPS, pa in PARAMS
+    dt = 0.01
+    tspan = TSPAN
+    t0, t1 = tspan
+    dynfun = dynode(iip, dp, dq)
+    p0, q0 = apa(iip, solution(t0, pa))
+    prob = DynamicalODEProblem(dynfun, p0, q0, tspan, pa)
+
+    if x || pa[1] == 1
+        sol = solve(prob, alg(); dt = dt)
+        calc = sol(t1)
+        # printerrors("$alg-$iip-$pa", calc, solution, pa, t1)
+        @test calc[1]≈solution(t1, pa)[1] rtol=errorbound(dt, d, calc[1])
+        @test calc[2]≈solution(t1, pa)[2] rtol=errorbound(dt, d, calc[2])
+    else
+        @test_throws ArgumentError solve(prob, alg(); dt = dt)
+    end
+end
+
+function motionfuncDirect1(dv,v,u,p,t)
+    # 1:Electron, 2: Be
+    ω_1,ω_2,γ,m_1,m_2,η,ω_d=p
+    dv[1]=-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1
+    dv[2]=-ω_2^2*u[2]-γ*u[1]/m_2
+end
+
+function motionfuncDirect1(v,u,p,t)
+    # 1:Electron, 2: Be
+    ω_1,ω_2,γ,m_1,m_2,η,ω_d=p
+    [-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1,-ω_2^2*u[2]-γ*u[1]/m_2]
+end
+
+param=[90386.15717208837, 3938.9288690708827, 8560.718748264337, 0.000544617021484666, 8.947079933513658, 0.7596480420227258, 78778.57738141765]
+u0_direct=zeros(2) # mm, mm
+v0_direct = [0.0, 135.83668926684385]
+tspan=(0.0, 1.321179076090661)
+prob_direct=SecondOrderODEProblem(motionfuncDirect1,v0_direct,u0_direct,tspan,param)
+dt=2e-8
+ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01)
+
+@testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS
+    sol=solve(prob_direct,alg(),dt=dt,saveat=0.01)
+    if alg <: Yoshida6
+        @test maximum(ref[4,:]-sol[4,:]) < 9e-3
+    else
+        @test maximum(ref[4,:]-sol[4,:]) < 3e-3
+    end
+end
diff --git a/test/downstream/delaydiffeq.jl b/test/downstream/delaydiffeq.jl
index 413fe86740..1c6d845c85 100644
--- a/test/downstream/delaydiffeq.jl
+++ b/test/downstream/delaydiffeq.jl
@@ -24,7 +24,7 @@ using Test
         @test sol.errors[:l∞] < error
 
         sol_scalar = solve(prob_scalar, ddealg)
-        @test sol.t ≈ sol_scalar.t
+        @test sol.t ≈ sol_scalar.t atol=1e-6
         @test sol[1, :] ≈ sol_scalar.u
     end
 end
diff --git a/test/gpu/Project.toml b/test/gpu/Project.toml
index 5a49c5b551..39a8540125 100644
--- a/test/gpu/Project.toml
+++ b/test/gpu/Project.toml
@@ -3,4 +3,4 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 
 [compat]
-CUDA = "4.2"
+CUDA = "5"
diff --git a/test/integrators/check_error.jl b/test/integrators/check_error.jl
index 9287e35dad..50c114e41c 100644
--- a/test/integrators/check_error.jl
+++ b/test/integrators/check_error.jl
@@ -5,7 +5,7 @@ u0 = 0.0 # explosion time is 1.0
 tspan = (0.0, 10.0)
 prob = ODEProblem(f_ec, u0, tspan)
 options = [:reltol => 1e-8, :abstol => 1e-8, :verbose => false]
-desired_code = ReturnCode.DtLessThanMin
+desired_code = ReturnCode.MaxIters
 
 # Test that sol.retcode is set to the correct value by various ways to
 # invoke integrator.
diff --git a/test/integrators/dae_initialization_tests.jl b/test/integrators/dae_initialization_tests.jl
index 77e13227bd..ae81100e0f 100644
--- a/test/integrators/dae_initialization_tests.jl
+++ b/test/integrators/dae_initialization_tests.jl
@@ -63,6 +63,37 @@ for alg in [Rodas5(autodiff = false), Trapezoid()]
     @test sum(sol[1]) ≈ 1
 end
 
+function rober_no_p(du, u, p, t)
+    y₁, y₂, y₃ = u
+    ;(k₁, k₂, k₃) = (0.04, 3e7, 1e4)
+    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
+    du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
+    du[3] = y₁ + y₂ + y₃ - 1
+    nothing
+end
+
+function rober_oop_no_p(du, u, p, t)
+    y₁, y₂, y₃ = u
+    ;(k₁, k₂, k₃) = (0.04, 3e7, 1e4)
+    du1 = -k₁ * y₁ + k₃ * y₂ * y₃
+    du2 = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
+    du3 = y₁ + y₂ + y₃ - 1
+    [du1, du2, du3]
+end
+
+# test oop and iip ODE initialization with parameters without eltype/length
+struct UnusedParam
+end
+for f in (ODEFunction(rober_no_p, mass_matrix = M), ODEFunction(rober_oop_no_p, mass_matrix = M))
+    local prob, probp
+    prob = ODEProblem(f, [1.0, 0.0, 1.0], (0.0, 1e5))
+    probp = ODEProblem(f, [1.0, 0.0, 1.0], (0.0, 1e5), UnusedParam)
+    for initializealg in (ShampineCollocationInit(), BrownFullBasicInit())
+        isapprox(init(prob, Rodas5(), abstol=1e-10; initializealg).u,
+                 init(prob, Rodas5(), abstol=1e-10; initializealg).u)
+    end
+end
+
 ## DAEProblem
 
 f = function (du, u, p, t)
@@ -170,6 +201,12 @@ integrator = init(prob, DABDF2(); initializealg = ShampineCollocationInit())
 @test_broken integrator.du[2]≈-1.0 atol=1e-9
 @test_broken integrator.u[3]≈2.0 atol=1e-9
 
+# test iip dae initialization with parameters without eltype/length
+probp = DAEProblem(f, du₀, u₀, tspan, UnusedParam(), differential_vars = differential_vars)
+for initializealg in (ShampineCollocationInit(), BrownFullBasicInit())
+    @test isapprox(init(probp, DABDF2(); initializealg).u, init(prob, DABDF2(); initializealg).u)
+end
+
 f = function (du, u, p, t)
     du - u
 end
@@ -196,6 +233,11 @@ integrator = init(prob, DABDF2())
 
 @test integrator.du[1]≈1.0 atol=1e-9
 @test integrator.du[2]≈1.0 atol=1e-9
+# test oop DAE initialization with parameters without eltype/length
+probp = DAEProblem(f, du₀, u₀, tspan, UnusedParam(), differential_vars = differential_vars)
+for initializealg in (ShampineCollocationInit(), BrownFullBasicInit())
+    @test isapprox(init(probp, DABDF2(); initializealg).u, init(prob, DABDF2(); initializealg).u)
+end
 
 # to test that we get the right NL solve we need a broken solver.
 struct BrokenNLSolve <: SciMLBase.AbstractNonlinearAlgorithm
@@ -216,3 +258,4 @@ prob = ODEProblem(f, ones(3), (0.0, 1.0))
 integrator = init(prob, Rodas5P(),
     initializealg = ShampineCollocationInit(1.0, BrokenNLSolve()))
 @test all(isequal(reinterpret(Float64, 0xDEADBEEFDEADBEEF)), integrator.u)
+
diff --git a/test/interface/ad_tests.jl b/test/interface/ad_tests.jl
index a8cd6ec1e6..1a85c97ad8 100644
--- a/test/interface/ad_tests.jl
+++ b/test/interface/ad_tests.jl
@@ -271,24 +271,6 @@ end
     @test ForwardDiff.derivative(g, 0.0)≈_u0 / 2 * exp(-0.5) rtol=rtol
 end
 
-function f(out, du, u, p, t)
-    out[1] = -p[1] * u[1] + 1e4 * u[2] * u[3] - du[1]
-    out[2] = +p[1] * u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
-    out[3] = u[1] + u[2] + u[3] - 1.0
-end
-
-p = [0.5]
-u₀ = [1.0, 0, 0]
-du₀ = [-0.04, 0.04, 0.0]
-tspan = (0.0, 100000.0)
-differential_vars = [true, true, false]
-prob = DAEProblem(f, du₀, u₀, tspan, p, differential_vars = differential_vars)
-
-function f(p)
-    sum(solve(remake(prob, p = p), DABDF2(), saveat = 0.1, abstol = 1e-6, reltol = 1e-6))
-end
-@test ForwardDiff.gradient(f, [0.5])[1]≈0 atol=1e-2
-
 # https://github.com/SciML/DifferentialEquations.jl/issues/903
 
 #ode function
diff --git a/test/interface/autodiff_error_tests.jl b/test/interface/autodiff_error_tests.jl
index c87589a428..e374166b54 100644
--- a/test/interface/autodiff_error_tests.jl
+++ b/test/interface/autodiff_error_tests.jl
@@ -34,12 +34,13 @@ tspan = (0.0, 1.0)
 prob = ODEProblem(lorenz!, u0, tspan)
 @test_throws OrdinaryDiffEq.FirstAutodiffJacError solve(prob, Rosenbrock23())
 
-function lorenz!(du, u, p, t)
+function lorenz2!(du, u, p, t)
     du[1] = 10.0(u[2] - u[1])
     a[1] = t
     du[2] = u[1] * (28.0 - u[3]) - u[2]
     du[3] = u[1] * u[2] - (8 / 3) * u[3]
 end
+prob = ODEProblem(lorenz2!, u0, tspan)
 @test_throws OrdinaryDiffEq.FirstAutodiffTgradError solve(prob, Rosenbrock23())
 
 ## Test that nothing is using duals when autodiff=false
diff --git a/test/interface/dae_ad_tests.jl b/test/interface/dae_ad_tests.jl
new file mode 100644
index 0000000000..d8976af291
--- /dev/null
+++ b/test/interface/dae_ad_tests.jl
@@ -0,0 +1,56 @@
+using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Test
+
+function rober(du, u, p, t)
+    y₁, y₂, y₃ = u
+    k₁, k₂, k₃ = p
+    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
+    du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
+    du[3] = y₁ + y₂ + y₃ - 1
+    nothing
+end
+function rober(u, p, t)
+    y₁, y₂, y₃ = u
+    k₁, k₂, k₃ = p
+    [-k₁ * y₁ + k₃ * y₂ * y₃,
+      k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2,
+      y₁ + y₂ + y₃ - 1]
+end
+M = [1.0 0 0
+    0 1.0 0
+    0 0 0]
+roberf = ODEFunction(rober, mass_matrix = M)
+roberf_oop = ODEFunction{false}(rober, mass_matrix = M)
+prob_mm = ODEProblem(roberf, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4))
+prob_mm_oop = ODEProblem(roberf_oop, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4))
+sol = solve(prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8)
+sol = solve(prob_mm_oop, Rodas5P(), reltol = 1e-8, abstol = 1e-8)
+
+function f(out, du, u, p, t)
+    out[1] = -p[1] * u[1] + p[3] * u[2] * u[3] - du[1]
+    out[2] = +p[1] * u[1] - p[2] * u[2]^2 - p[3] * u[2] * u[3] - du[2]
+    out[3] = u[1] + u[2] + u[3] - 1.0
+end
+function f(du, u, p, t)
+    [-p[1] * u[1] + p[3] * u[2] * u[3] - du[1],
+     +p[1] * u[1] - p[2] * u[2]^2 - p[3] * u[2] * u[3] - du[2],
+     u[1] + u[2] + u[3] - 1.0]
+end
+p = [0.04, 3e7, 1e4]
+u₀ = [1.0, 0, 0]
+du₀ = [-0.04, 0.04, 0.0]
+tspan = (0.0, 100000.0)
+differential_vars = [true, true, false]
+prob = DAEProblem(f, du₀, u₀, tspan, p, differential_vars = differential_vars)
+prob_oop = DAEProblem{false}(f, du₀, u₀, tspan, p, differential_vars = differential_vars)
+sol1 = solve(prob, DFBDF(), dt=1e-5, abstol = 1e-8, reltol = 1e-8)
+
+# These tests flex differentiation of the solver and through the initialization
+# To only test the solver part and isolate potential issues, set the initialization to consistent
+@testset "Inplace: $(isinplace(_prob)), DAEProblem: $(_prob isa DAEProblem), BrownBasic: $(initalg isa BrownFullBasicInit), Autodiff: $autodiff" for _prob in [prob, prob_mm, prob_oop, prob_mm_oop], initalg in [BrownFullBasicInit(), ShampineCollocationInit()], autodiff in [true, false]
+    alg = _prob isa DAEProblem ? DFBDF(;autodiff) : Rodas5P(;autodiff)
+    function f(p)
+        sol = solve(remake(_prob, p = p), alg, abstol = 1e-14, reltol = 1e-14, initializealg = initalg)
+        sum(sol)
+    end
+    @test ForwardDiff.gradient(f, [0.04, 3e7, 1e4]) ≈ [0,0,0] atol=1e-8
+end
diff --git a/test/interface/interpolation_output_types.jl b/test/interface/interpolation_output_types.jl
new file mode 100644
index 0000000000..f16a0fcdee
--- /dev/null
+++ b/test/interface/interpolation_output_types.jl
@@ -0,0 +1,41 @@
+using OrdinaryDiffEq, RecursiveArrayTools, Test
+
+# in terms of the voltage across all three elements
+rlc1!(v′,v,(R,L,C),t) = -(v′/R + v/L)/C
+identity_f(v,u,p,t) = v # needed to form second order dynamical ODE
+
+setup_rlc(R,L,C;v_init=0.0,v′_init=0.0,tspan=(0.0,50.0)) =
+    DynamicalODEProblem{false}(rlc1!,identity_f,v′_init,v_init,tspan,(R,L,C))
+
+# simulate voltage impulse
+R,L,C = 10, 0.3, 2
+
+prob = setup_rlc(R,L,C,v_init=2.0)
+
+res1 = solve(prob,Vern8(),dt=1/10,saveat=1/10)
+res3 = solve(prob,CalvoSanz4(),dt=1/10,saveat=1/10)
+
+sol = solve(prob,CalvoSanz4(),dt=1/10)
+@test sol(0.32) isa OrdinaryDiffEq.ArrayPartition
+@test sol(0.32, Val{1}) isa OrdinaryDiffEq.ArrayPartition
+@test sol(0.32, Val{2}) isa OrdinaryDiffEq.ArrayPartition
+@test sol(0.32, Val{3}) isa OrdinaryDiffEq.ArrayPartition
+
+function f(du,u,p,t)
+    du .= u
+    nothing
+end
+dprob = DiscreteProblem(f, [1,2,3], (0,100))
+sol = solve(dprob, FunctionMap())
+@test sol(0:0.1:100;idxs=1) isa DiffEqArray
+@test length(sol(0:0.1:100;idxs=1)) == length(0:0.1:100)
+@test length(sol(0:0.1:100;idxs=1).u[1]) == 1
+sol(0:0.1:100;idxs=[1,2])
+
+@test sol(0:0.1:100;idxs=[1,2]) isa DiffEqArray
+@test length(sol(0:0.1:100;idxs=[1,2])) == length(0:0.1:100)
+@test length(sol(0:0.1:100;idxs=[1,2]).u[1]) == 2
+
+@test sol(0:0.1:100) isa DiffEqArray
+@test length(sol(0:0.1:100)) == length(0:0.1:100)
+@test length(sol(0:0.1:100).u[1]) == 3
diff --git a/test/interface/ode_initdt_tests.jl b/test/interface/ode_initdt_tests.jl
index 04f8bd8f04..9e75bbf252 100644
--- a/test/interface/ode_initdt_tests.jl
+++ b/test/interface/ode_initdt_tests.jl
@@ -28,7 +28,7 @@ prob = remake(prob, u0 = u0, tspan = tspan)
 tspan = T.((2000, 2100))
 prob = remake(prob, tspan = tspan)
 # set maxiters to prevent infinite loop on test failure
-@test_throws ArgumentError solve(prob, Euler(); dt = T(0.0001), maxiters = 10)
+@test solve(prob, Euler(); dt = T(0.0001), maxiters = 10).retcode == SciMLBase.ReturnCode.MaxIters
 
 function rober(du, u, p, t)
     y₁, y₂, y₃ = u
diff --git a/test/interface/ode_saveat_tests.jl b/test/interface/ode_saveat_tests.jl
index 858fada876..5e1a6de97e 100644
--- a/test/interface/ode_saveat_tests.jl
+++ b/test/interface/ode_saveat_tests.jl
@@ -160,7 +160,7 @@ integ = init(ODEProblem((u, p, t) -> u, 0.0, (0.0, 1.0)), Tsit5(), saveat = _sav
     save_end = false)
 add_tstop!(integ, 2.0)
 solve!(integ)
-@test integ.sol.t == _saveat
+@test integ.sol.t == _saveat[1:end-1]
 
 # Catch save for maxiters
 ode = ODEProblem((u, p, t) -> u, 1.0, (0.0, 1.0))
@@ -187,3 +187,28 @@ prob = ODEProblem(SIR!, [0.99, 0.01, 0.0], (t_obs[1], t_obs[end]), [0.20, 0.15])
 sol = solve(prob, DP5(), reltol = 1e-6, abstol = 1e-6, saveat = t_obs)
 @test maximum(sol) <= 1
 @test minimum(sol) >= 0
+
+@testset "Proper save_start and save_end behavior" begin
+    function f2(du, u, p, t)
+        du[1] = -cos(u[1]) * u[1]
+    end
+    prob = ODEProblem(f2, [10], (0.0, 0.4))
+
+    @test solve(prob, Tsit5(); saveat = 0:.1:.4).t == [0.0; 0.1; 0.2; 0.3; 0.4]
+    @test solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = true, save_end = true).t == [0.0; 0.1; 0.2; 0.3; 0.4]
+    @test solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = false, save_end = false).t == [0.1; 0.2; 0.3]
+
+    ts = solve(prob, Tsit5()).t
+    @test 0.0 in ts
+    @test 0.4 in ts
+    ts = solve(prob, Tsit5(); save_start = true, save_end = true).t
+    @test 0.0 in ts
+    @test 0.4 in ts
+    ts = solve(prob, Tsit5(); save_start = false, save_end = false).t
+    @test 0.0 ∉ ts
+    @test 0.4 ∉ ts
+
+    @test solve(prob, Tsit5(); saveat = [.2]).t == [0.2]
+    @test solve(prob, Tsit5(); saveat = [.2], save_start = true, save_end = true).t == [0.0; 0.2; 0.4]
+    @test solve(prob, Tsit5(); saveat = [.2], save_start = false, save_end = false).t == [0.2]
+end
diff --git a/test/interface/static_array_tests.jl b/test/interface/static_array_tests.jl
index cae40bc4f4..ff2bcbb3cf 100644
--- a/test/interface/static_array_tests.jl
+++ b/test/interface/static_array_tests.jl
@@ -2,7 +2,7 @@ using StaticArrays, Test
 using OrdinaryDiffEq
 using RecursiveArrayTools
 
-u0 = [fill(2, MVector{2, Float64}), ones(MVector{2, Float64})]
+u0 = VectorOfArray([fill(2, MVector{2, Float64}), ones(MVector{2, Float64})])
 g(u, p, t) = SA[u[1] + u[2], u[1]]
 f = (du, u, p, t) -> begin
     for i in 1:2
@@ -17,7 +17,7 @@ sol = solve(ode, Tsit5())
 sol = solve(ode, Vern9())
 @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
 
-u0 = [fill(2, SVector{2, Float64}), ones(SVector{2, Float64})]
+u0 = VectorOfArray([fill(2, SVector{2, Float64}), ones(SVector{2, Float64})])
 ode = ODEProblem(f, u0, (0.0, 1.0))
 sol = solve(ode, Euler(), dt = 1e-2)
 @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
diff --git a/test/regression/inference.jl b/test/regression/inference.jl
index dc28b6bc33..9e62015f54 100644
--- a/test/regression/inference.jl
+++ b/test/regression/inference.jl
@@ -7,13 +7,13 @@ using Test
 
     inferred = [BS3(), Tsit5(), RK4(), Vern6()]
     for alg in inferred
-        @test_broken init(prob, alg)
-        @test_broken init(prob2D, alg)
+        @inferred init(prob, alg)
+        @inferred init(prob2D, alg)
     end
 
-    notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()]
-    for alg in notinferred
-        @test_broken @inferred init(prob, alg)
-        @test_broken @inferred init(prob2D, alg)
-    end
+    #notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()]
+    #for alg in notinferred
+    #    @test_broken @inferred init(prob, alg)
+    #    @test_broken @inferred init(prob2D, alg)
+    #end
 end
diff --git a/test/regression/ode_dense_tests.jl b/test/regression/ode_dense_tests.jl
index dc60c17c0b..5eb1441d6b 100644
--- a/test/regression/ode_dense_tests.jl
+++ b/test/regression/ode_dense_tests.jl
@@ -537,10 +537,10 @@ regression_test(Rodas4P(), 4e-5, 6e-5, test_diff1 = true, nth_der = 1, dertol =
 regression_test(Rodas4P2(), 2e-5, 3e-5, test_diff1 = true, nth_der = 1, dertol = 1e-13)
 
 # Rodas5
-regression_test(Rodas5(), 2e-6, 3e-6, test_diff1 = true, nth_der = 1, dertol = 1e-13)
+regression_test(Rodas5(), 2e-6, 3e-6, test_diff1 = true, nth_der = 3, dertol = 5e-1)
 
 # Rodas5P
-regression_test(Rodas5P(), 2e-5, 3e-5, test_diff1 = true, nth_der = 1, dertol = 1e-13)
+regression_test(Rodas5P(), 2e-5, 3e-5, test_diff1 = true, nth_der = 3, dertol = 5e-1)
 
 # ExplicitRK
 regression_test(ExplicitRK(), 7e-5, 2e-4)
diff --git a/test/runtests.jl b/test/runtests.jl
index 395c139a84..183e3a4b6f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -44,6 +44,7 @@ end
         @time @safetestset "Complex Tests" include("interface/complex_tests.jl")
         @time @safetestset "Ndim Complex Tests" include("interface/ode_ndim_complex_tests.jl")
         @time @safetestset "Number Type Tests" include("interface/ode_numbertype_tests.jl")
+        @time @safetestset "Interpolation Output Type Tests" include("interface/interpolation_output_types.jl")
         @time @safetestset "Stiffness Detection Tests" include("interface/stiffness_detection_test.jl")
         @time @safetestset "Composite Interpolation Tests" include("interface/composite_interpolation.jl")
         @time @safetestset "Export tests" include("interface/export_tests.jl")
@@ -83,8 +84,9 @@ end
     end
 
     if !is_APPVEYOR && (GROUP == "All" || GROUP == "InterfaceV" || GROUP == "Interface")
-        @time @safetestset "AD Tests" include("interface/interpolation_derivative_error_tests.jl")
+        @time @safetestset "Interpolation Derivative Error Tests" include("interface/interpolation_derivative_error_tests.jl")
         @time @safetestset "AD Tests" include("interface/ad_tests.jl")
+        @time @safetestset "DAE AD Tests" include("interface/dae_ad_tests.jl")
         @time @safetestset "Newton Tests" include("interface/newton_tests.jl")
         @time @safetestset "DAE Initialize Integration" include("interface/dae_initialize_integration.jl")
     end