diff --git a/Project.toml b/Project.toml
index 868f2451..ebe9c3b2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -48,7 +48,6 @@ QuadGK = "2.9"
 Reexport = "1.0"
 SafeTestsets = "0.1"
 SciMLBase = "2.6"
-SciMLSensitivity = "7.41"
 StaticArrays = "1"
 Test = "1"
 Zygote = "0.6.60"
@@ -66,10 +65,9 @@ FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Arblib", "SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"]
+test = ["Aqua", "Arblib", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"]
diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl
index 5aec07c3..300b6063 100644
--- a/ext/IntegralsForwardDiffExt.jl
+++ b/ext/IntegralsForwardDiffExt.jl
@@ -3,7 +3,8 @@ using Integrals
 isdefined(Base, :get_extension) ? (using ForwardDiff) : (using ..ForwardDiff)
 ### Forward-Mode AD Intercepts
 
-# Direct AD on solvers with QuadGK and HCubature
+#= Direct AD on solvers with QuadGK and HCubature
+# incompatible with iip since types must change
 function Integrals.__solvebp(cache, alg::QuadGKJL, sensealg, domain,
         p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N};
         kwargs...) where {T, V, P, N}
@@ -15,11 +16,14 @@ function Integrals.__solvebp(cache, alg::HCubatureJL, sensealg, domain,
         kwargs...) where {T, V, P, N}
     Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...)
 end
+=#
+
+# TODO: add the pushforward for derivative w.r.t lb, and ub (and then combinations?)
 
 # Manually split for the pushforward
 function Integrals.__solvebp(cache, alg, sensealg, domain,
-        p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N};
-        kwargs...) where {T, V, P, N}
+        p::Union{D,AbstractArray{<:D}};
+        kwargs...) where {T, V, P, D<:ForwardDiff.Dual{T, V, P}}
 
     # we need the output type to avoid perturbation confusion while unwrapping nested duals
     # We compute a vector-valued integral of the primal and dual simultaneously
@@ -32,11 +36,11 @@ function Integrals.__solvebp(cache, alg, sensealg, domain,
             len,
             size(cache.f.integrand_prototype)...)
 
-        dfdp_ = function (out, x, p)
-            dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p)
+        dfdp_ = function (out, x, _p)
+            dualp = reinterpret(ForwardDiff.Dual{T, V, P}, _p)
             dout = reinterpret(reshape, DT, out)
-            cache.f(dout, x, dualp)
-            return out
+            cache.f(dout, x, p isa D ? only(dualp) : reshape(dualp, size(p)))
+            return
         end
         dfdp = if cache.f isa BatchIntegralFunction
             BatchIntegralFunction{true}(dfdp_, dual_prototype)
@@ -55,9 +59,9 @@ function Integrals.__solvebp(cache, alg, sensealg, domain,
         DT = y isa AbstractArray ? eltype(y) : typeof(y)
         elt = unwrap_dualvaltype(DT)
 
-        dfdp_ = function (x, p)
-            dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p)
-            ys = cache.f(x, dualp)
+        dfdp_ = function (x, _p)
+            dualp = reinterpret(ForwardDiff.Dual{T, V, P}, _p)
+            ys = cache.f(x, p isa D ? only(dualp) : reshape(dualp, size(p)))
             ys_ = ys isa AbstractArray ? ys : [ys]
             # we need to reshape in order for batching to be consistent
             return reinterpret(reshape, elt, ys_)
@@ -70,7 +74,7 @@ function Integrals.__solvebp(cache, alg, sensealg, domain,
     end
 
     ForwardDiff.can_dual(elt) || ForwardDiff.throw_cannot_dual(elt)
-    rawp = copy(reinterpret(V, p))
+    rawp = p isa D ? reinterpret(V, [p]) : copy(reinterpret(V, vec(p)))
 
     prob = Integrals.build_problem(cache)
     dp_prob = remake(prob, f = dfdp, p = rawp)
diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl
index d5dcc034..90d6f882 100644
--- a/ext/IntegralsZygoteExt.jl
+++ b/ext/IntegralsZygoteExt.jl
@@ -1,4 +1,5 @@
 module IntegralsZygoteExt
+using LinearAlgebra: dot
 using Integrals
 if isdefined(Base, :get_extension)
     using Zygote
@@ -11,6 +12,7 @@ else
 end
 ChainRulesCore.@non_differentiable Integrals.checkkwargs(kwargs...)
 ChainRulesCore.@non_differentiable Integrals.isinplace(f, args...)    # fixes #99
+ChainRulesCore.@non_differentiable Integrals.init_cacheval(alg, prob)
 
 function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, sensealg, domain,
         p;
@@ -24,23 +26,25 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
         # https://juliadiff.org/ChainRulesCore.jl/dev/design/many_tangents.html#manytypes
         if isinplace(cache)
             # zygote doesn't support mutation, so we build an oop pullback
-            dx = similar(cache.f.integrand_prototype)
-            _f = x -> cache.f(dx, x, p)
             if sensealg.vjp isa Integrals.ZygoteVJP
                 if cache.f isa BatchIntegralFunction
+                    dx = similar(cache.f.integrand_prototype, size(cache.f.integrand_prototype)[begin:end-1]..., 1)
+                    _f = x -> (cache.f(dx, x, p); dx)
                     # TODO: let the user pass a batched jacobian so we can return a BatchIntegralFunction
                     dfdp_ = function (x, p)
                         x_ = x isa AbstractArray ? reshape(x, size(x)..., 1) : [x]
                         z, back = Zygote.pullback(p) do p
-                            _dx = Zygote.Buffer(dx, size(dx)[begin:(end - 1)]..., 1)
+                            _dx = Zygote.Buffer(dx)
                             cache.f(_dx, x_, p)
                             copy(_dx)
                         end
                         return back(z .= (Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) :
-                                          [Δ]))[1]
+                                          Δ))[1]
                     end
                     dfdp = IntegralFunction{false}(dfdp_, nothing)
                 else
+                    dx = similar(cache.f.integrand_prototype)
+                    _f = x -> (cache.f(dx, x, p); dx)
                     dfdp_ = function (x, p)
                         _, back = Zygote.pullback(p) do p
                             _dx = Zygote.Buffer(dx)
@@ -62,8 +66,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
                     dfdp_ = function (x, p)
                         x_ = x isa AbstractArray ? reshape(x, size(x)..., 1) : [x]
                         z, back = Zygote.pullback(p -> cache.f(x_, p), p)
-                        return back(z .= (Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) :
-                                          [Δ]))[1]
+                        return back(Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) : [Δ])[1]
                     end
                     dfdp = IntegralFunction{false}(dfdp_, nothing)
                 else
@@ -98,13 +101,15 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
 
         lb, ub = domain
         if lb isa Number
-            dlb = cache.f isa BatchIntegralFunction ? -_f([lb]) : -_f(lb)
-            dub = cache.f isa BatchIntegralFunction ? _f([ub]) : _f(ub)
+            # TODO replace evaluation at endpoint (which anyone can do without Integrals.jl)
+            # with integration of dfdx uing the same quadrature
+            dlb = cache.f isa BatchIntegralFunction ? -batch_unwrap(_f([lb])) : -_f(lb)
+            dub = cache.f isa BatchIntegralFunction ?  batch_unwrap(_f([ub])) : _f(ub)
             return (NoTangent(),
                 NoTangent(),
                 NoTangent(),
                 NoTangent(),
-                Tangent{typeof(domain)}(dlb, dub),
+                Tangent{typeof(domain)}(dot(dlb, Δ), dot(dub, Δ)),
                 dp)
         else
             # we need to compute 2*length(lb) integrals on the faces of the hypercube, as we
@@ -123,6 +128,8 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
     out, quadrature_adjoint
 end
 
+batch_unwrap(x::AbstractArray) = dropdims(x; dims=ndims(x))
+
 Zygote.@adjoint function Zygote.literal_getproperty(sol::SciMLBase.IntegralSolution,
         ::Val{:u})
     sol.u, Δ -> (SciMLBase.build_solution(sol.prob, sol.alg, Δ, sol.resid),)
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index 87b326f7..061549d8 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -1,243 +1,309 @@
 using Integrals, Zygote, FiniteDiff, ForwardDiff#, SciMLSensitivity
 using Cuba, Cubature
+using FastGaussQuadrature
 using Test
 
-### One Dimensional
-
-f(x, p) = sum(sin.(p[1] * x))
-lb = 1.0
-ub = 3.0
-p = 2.0
-prob = IntegralProblem(f, lb, ub, p)
-sol = solve(prob, QuadGKJL(), reltol = 1e-3, abstol = 1e-3)
-
-function testf(lb, ub, p)
-    prob = IntegralProblem(f, lb, ub, p)
-    sin(solve(prob, QuadGKJL(), reltol = 1e-3, abstol = 1e-3)[1])
+max_dim_test = 2
+max_nout_test = 2
+
+reltol = 1e-3
+abstol = 1e-3
+
+alg_req = Dict(
+    QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
+        allows_iip = false),
+    QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1,
+        allows_iip = true),
+    HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1,
+        max_dim = Inf, allows_iip = true),
+    # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf,
+        # allows_iip = true),
+    CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1,
+        max_dim = Inf, allows_iip = true),
+    CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1,
+        max_dim = Inf, allows_iip = true),
+    # CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
+    #     allows_iip = true),
+    # CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
+    #     allows_iip = true),
+    # CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2,
+    #     max_dim = Inf, allows_iip = true),
+    # CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf,
+    #     allows_iip = true),
+)
+
+# integrands should have same shape as parameters, independent of dimensionality
+integrands =(
+    (x, p) -> map(q -> prod(y -> sin(y*q), x), p),
+)
+
+# function to turn the output into a scalar / test different tangent types
+scalarize_solution = (
+    sol -> sin(sum(sol)),
+    sol -> sin(sol[1]),
+)
+
+# we will be able to use broadcasting for this after https://github.com/FluxML/Zygote.jl/pull/1488
+function buffer_copyto!(y, x)
+    for (j, i) in zip(eachindex(y), eachindex(x))
+        y[j] = x[i]
+    end
+    return y
+end
+function f_helper!(f, y, x, p)
+    buffer_copyto!(y, f(x, p))
+    return
 end
-dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p)
-dlb2 = FiniteDiff.finite_difference_derivative(lb -> testf(lb, ub, p), lb)
-dub2 = FiniteDiff.finite_difference_derivative(ub -> testf(lb, ub, p), ub)
-dp2 = FiniteDiff.finite_difference_derivative(p -> testf(lb, ub, p), p)
-
-# dlb3 = ForwardDiff.derivative(lb->testf(lb,ub,p),lb)
-# dub3 = ForwardDiff.derivative(ub->testf(lb,ub,p),ub)
-dp3 = ForwardDiff.derivative(p -> testf(lb, ub, p), p)
-
-#@test dlb1 ≈ dlb3
-#@test dub1 ≈ dub3
-@test dp1 ≈ dp3
-
-### N-dimensional
-
-f(x, p) = sum(sin.(x .* p))
-lb = ones(2)
-ub = 3ones(2)
-p = [1.5, 2.0]
-prob = IntegralProblem(f, lb, ub, p)
-sol = solve(prob, CubaCuhre(), reltol = 1e-3, abstol = 1e-3)
 
-function testf(lb, ub, p)
-    prob = IntegralProblem(f, lb, ub, p)
-    sin(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1])
+# the Zygote implementation is inconsistent about 0-d so we hijack it
+struct Scalar{T<:Real} <: Real
+    x::T
+end
+Base.iterate(a::Scalar) = (a.x, nothing)
+Base.iterate(::Scalar, _) = nothing
+Base.IteratorSize(::Type{Scalar{T}}) where {T} = Base.HasShape{0}()
+Base.eltype(::Type{Scalar{T}}) where {T} = T
+Base.length(a::Scalar) = 1
+Base.size(::Scalar) = ()
+Base.:+(a::Scalar, b::Scalar) = Scalar(a.x+b.x)
+Base.:*(a::Number, b::Scalar) = a*b.x
+Base.:*(a::Scalar, b::Number) = a.x*b
+Base.:*(a::Scalar, b::Scalar) = Scalar(a.x*b.x)
+Base.zero(a::Scalar) = Scalar(zero(a.x))
+Base.map(f, a::Scalar) = map(f, a.x)
+(::Type{T})(a::Scalar) where {T<:Real} = T(a.x)
+struct ScalarAxes end # the implementation doesn't preserve singleton axes
+Base.axes(::Scalar) = ScalarAxes()
+Base.iterate(::ScalarAxes) = nothing
+Base.reshape(A::AbstractArray, ::ScalarAxes) = Scalar(only(A))
+
+# here we assume f evaluated at scalar inputs gives a scalar output
+# p will be able to be a number  after https://github.com/FluxML/Zygote.jl/pull/1489
+# p will be able to be a 0-array after https://github.com/FluxML/Zygote.jl/pull/1491
+# p can't be either without both prs
+function batch_helper(f, x, p)
+    t = f(zero(eltype(x)), zero(eltype(eltype(p))))
+    typeof(t).([f(y, q) for q in p, y in eachslice(x; dims=ndims(x))])
 end
 
-function testf2(lb, ub, p)
-    prob = IntegralProblem(f, lb, ub, p)
-    sin(solve(prob, HCubatureJL(), reltol = 1e-6, abstol = 1e-6)[1])
+function batch_helper!(f, y, x, p)
+    buffer_copyto!(y, batch_helper(f, x, p))
+    return
 end
 
-dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p)
-dlb2 = FiniteDiff.finite_difference_gradient(lb -> testf(lb, ub, p), lb)
-dub2 = FiniteDiff.finite_difference_gradient(ub -> testf(lb, ub, p), ub)
-dp2 = FiniteDiff.finite_difference_gradient(p -> testf(lb, ub, p), p)
+# helper function / test runner
+do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol)
+    testf = function (lb, ub, p)
+        prob = IntegralProblem(f, (lb, ub), p)
+        scalarize(solve(prob, alg; reltol, abstol))
+    end
+    testf(lb, ub, p)
 
-@test_broken dlb1 ≈ dlb2
-@test_broken dub1 ≈ dub2
-@test dp1 ≈ dp2
+    dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p isa Number && f isa BatchIntegralFunction ? Scalar(p) : p)
 
-# dlb3 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
-# dub3 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
-dp3 = ForwardDiff.gradient(p -> testf2(lb, ub, p), p)
+    f_lb = lb -> testf(lb, ub, p)
+    f_ub = ub -> testf(lb, ub, p)
 
-#@test dlb1 ≈ dlb3
-#@test dub1 ≈ dub3
-@test dp1 ≈ dp3
+    dlb = lb isa AbstractArray ? :gradient : :derivative
+    dub = ub isa AbstractArray ? :gradient : :derivative
 
-# dlb4 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
-# dub4 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
-dp4 = ForwardDiff.gradient(p -> testf(lb, ub, p), p)
+    dlb2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dlb))(f_lb, lb)
+    dub2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dub))(f_ub, ub)
 
-#@test dlb1 ≈ dlb4
-#@test dub1 ≈ dub4
-@test dp1 ≈ dp4
+    if lb isa Number
+        @test dlb1 ≈ dlb2 atol=abstol rtol=reltol
+        @test dub1 ≈ dub2 atol=abstol rtol=reltol
+    else # TODO: implement multivariate limit derivatives in ZygoteExt
+        @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol
+        @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol
+    end
 
-### N-dimensional N-out
+    # TODO: implement limit derivatives in ForwardDiffExt
+    @test_broken dlb2 ≈ getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol
+    @test_broken dub2 ≈ getproperty(ForwardDiff, dub)(dfdub, ub) atol=abstol rtol=reltol
 
-f(x, p) = sin.(x .* p)
-lb = ones(2)
-ub = 3ones(2)
-p = [1.5, 2.0]
-prob = IntegralProblem(f, lb, ub, p, nout = 2)
-sol = solve(prob, CubaCuhre(), reltol = 1e-3, abstol = 1e-3)
+    f_p = p -> testf(lb, ub, p)
 
-function testf(lb, ub, p)
-    prob = IntegralProblem(f, lb, ub, p, nout = 2)
-    sin(sum(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)))
-end
+    dp = p isa AbstractArray ? :gradient : :derivative
 
-function testf2(lb, ub, p)
-    prob = IntegralProblem(f, lb, ub, p, nout = 2)
-    sin(sum(solve(prob, HCubatureJL(), reltol = 1e-6, abstol = 1e-6)))
-end
+    dp2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dp))(f_p, p)
+    dp3 = getproperty(ForwardDiff, dp)(f_p, p)
 
-dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p)
-dlb2 = FiniteDiff.finite_difference_gradient(lb -> testf(lb, ub, p), lb)
-dub2 = FiniteDiff.finite_difference_gradient(ub -> testf(lb, ub, p), ub)
-dp2 = FiniteDiff.finite_difference_gradient(p -> testf(lb, ub, p), p)
+    @test dp1 ≈ dp2 atol=abstol rtol=reltol
+    @test dp2 ≈ dp3 atol=abstol rtol=reltol
 
-@test_broken dlb1 ≈ dlb2
-@test_broken dub1 ≈ dub2
-@test dp1 ≈ dp2
-
-# dlb3 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
-# dub3 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
-dp3 = ForwardDiff.gradient(p -> testf2(lb, ub, p), p)
-
-#@test dlb1 ≈ dlb3
-#@test dub1 ≈ dub3
-@test dp1 ≈ dp3
+    return
+end
 
-# dlb4 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
-# dub4 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
-dp4 = ForwardDiff.gradient(p -> testf(lb, ub, p), p)
 
-#@test dlb1 ≈ dlb4
-#@test dub1 ≈ dub4
-@test dp1 ≈ dp4
+### One Dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
 
-### Batch Single dim
-f(x, p) = x * p[1] .+ p[2] * p[3]  # scalar integrand
+    @info "One-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i
+    do_tests(; f, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
+end
 
-lb = 1.0
-ub = 3.0
-p = [2.0, 3.0, 4.0]
-prob = IntegralProblem(f, lb, ub, p)
+## One-dimensional nout
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
 
-function testf3(lb, ub, p; f = f)
-    prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 1)
-    solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)[1]
+    @info "One-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i nout
+    do_tests(; f, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
-dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p)
-dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1] # TODO fix: LoadError: DimensionMismatch("variable with size(x) == (1, 15) cannot have a gradient with size(dx) == (15,)")
-dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p)
+### N-dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
 
-@test dp1 ≈ dp3 #passes
-@test dp2 ≈ dp3 #passes
+    @info "Multi-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i dim
+    do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
+end
 
-### Batch single dim, nout
-f(x, p) = (x' * p[1] .+ p[2] * p[3]) .* [1; 2]
+### N-dimensional nout
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
 
-lb = 1.0
-ub = 3.0
-p = [2.0, 3.0, 4.0]
-prob = IntegralProblem(f, lb, ub, p)
+    @info "Multi-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i dim nout
+    do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
+end
 
-function testf3(lb, ub, p; f = f)
-    prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 2)
-    sum(solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3))
+### One Dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
+
+    @info "One-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i
+    fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1))
+    do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
 end
 
-dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p)
-dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1]
-dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p)
+## One-dimensional nout
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
 
-@test dp1 ≈ dp3 #passes
-@test dp2 ≈ dp3 #passes
+    @info "One-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i nout
+    fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout))
+    do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
+end
 
-### Batch multi dim
-f(x, p) = x[1, :] * p[1] .+ p[2] * p[3]
+### N-dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
+
+    @info "Multi-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i dim
+    fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1))
+    do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
+end
 
-lb = [1.0, 1.0]
-ub = [3.0, 3.0]
-p = [2.0, 3.0, 4.0]
-prob = IntegralProblem(f, lb, ub, p)
+### N-dimensional nout iip
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
 
-function testf3(lb, ub, p; f = f)
-    prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 1)
-    solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)[1]
+    @info "Multi-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i dim nout
+    fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout))
+    do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
-dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p)
-dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1]
-dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p)
+### Batch, One Dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+    req.allows_batch || continue
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
 
-@test dp1 ≈ dp3
-@test dp2 ≈ dp3
-
-### Batch multi dim nout
-f(x, p) = x * p[1] .+ p[2] * p[3]
+    @info "Batched, one-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i
+    bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
+    do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
+end
 
-lb = [1.0, 1.0]
-ub = [3.0, 3.0]
-p = [2.0, 3.0, 4.0]
-prob = IntegralProblem(f, lb, ub, p)
+## Batch, One-dimensional nout
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+    req.allows_batch || continue
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
 
-function testf3(lb, ub, p; f = f)
-    prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 2)
-    sum(solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3))
+    @info "Batched, one-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i nout
+    bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
+    do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
-dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p)
-dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1]
-dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p)
+### Batch, N-dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+    req.allows_batch || continue
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
 
-@test dp1 ≈ dp3
-@test dp2 ≈ dp3
-
-## iip Batch multi dim
-function g(dx, x, p)
-    dx .= sum(x * p[1] .+ p[2] * p[3], dims = 1)
+    @info "Batched, multi-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i dim
+    bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
+    do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
 end
 
-lb = [1.0, 1.0]
-ub = [3.0, 3.0]
-p = [2.0, 3.0, 4.0]
-prob = IntegralProblem(g, lb, ub, p)
+### Batch, N-dimensional nout
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+    req.allows_batch || continue
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
 
-function testf3(lb, ub, p; f = g)
-    prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 1)
-    solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)[1]
+    @info "Batch, multi-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i dim nout
+    bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
+    do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
-testf3(lb, ub, p)
+### Batch, one-dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+    req.allows_batch || continue
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
 
-dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p)
-dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1]
-dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p)
+    @info "Batched, one-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i
+    bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0))
+    do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
+end
 
-@test dp1 ≈ dp3
-@test dp2 ≈ dp3
+## Batch, one-dimensional nout
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+    req.allows_batch || continue
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= 1 || continue
 
-## iip Batch mulit dim nout
-function g(dx, x, p)
-    dx .= x * p[1] .+ p[2] * p[3]
+    @info "Batched, one-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i nout
+    bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0))
+    do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
-lb = [1.0, 1.0]
-ub = [3.0, 3.0]
-p = [2.0, 3.0, 4.0]
-prob = IntegralProblem(g, lb, ub, p)
+### Batch, N-dimensional
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+    req.allows_batch || continue
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
 
-function testf3(lb, ub, p; f = g)
-    prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 2)
-    sum(solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3))
+    @info "Batched, multi-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i dim
+    bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0))
+    do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
 end
 
-dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p)
-dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1]
-dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p)
+### Batch, N-dimensional nout iip
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+    req.allows_batch || continue
+    req.allows_iip || continue
+    req.nout > 1 || continue
+    req.min_dim <= dim <= req.max_dim || continue
 
-@test dp1 ≈ dp3
-@test dp2 ≈ dp3
+    @info "Batched, multi-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i dim nout
+    bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0))
+    do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
+end