diff --git a/DifferentiationInterface/ext/DifferentiationInterfaceStaticArraysExt/DifferentiationInterfaceStaticArraysExt.jl b/DifferentiationInterface/ext/DifferentiationInterfaceStaticArraysExt/DifferentiationInterfaceStaticArraysExt.jl
index f450aaadc..eb8157d99 100644
--- a/DifferentiationInterface/ext/DifferentiationInterfaceStaticArraysExt/DifferentiationInterfaceStaticArraysExt.jl
+++ b/DifferentiationInterface/ext/DifferentiationInterfaceStaticArraysExt/DifferentiationInterfaceStaticArraysExt.jl
@@ -2,14 +2,52 @@ module DifferentiationInterfaceStaticArraysExt
 
 using ADTypes: AutoForwardDiff, AutoEnzyme
 import DifferentiationInterface as DI
-using StaticArrays: SArray, StaticArray
+using StaticArrays: SArray, Size, SOneTo, StaticArray, similar_type
+
+# a flat vector of the same kind (e.g. `MVector` for an `MMatrix`) keeps concatenated blocks `SArray`s or `MArray`s
+DI.vec_block(a::StaticArray) = similar_type(a, Size(length(a)))(a)
 
 function DI.stack_vec_col(t::NTuple{B, <:StaticArray}) where {B}
-    return hcat(map(vec, t)...)
+    return hcat(map(DI.vec_block, t)...)
 end
 
 function DI.stack_vec_row(t::NTuple{B, <:StaticArray}) where {B}
-    return vcat(transpose.(map(vec, t))...)
+    return vcat(transpose.(map(DI.vec_block, t))...)
+end
+
+# The number of seeds `N` and of batches `A` are deduced from the seed type,
+# so that concatenating the blocks preserves static sizes.
+
+function static_batch_count(::Type{T}, B::Int) where {T <: StaticArray}
+    N = length(T)
+    A = cld(N, B)
+    # otherwise `hcat()` would silently return an empty `Vector{Any}`
+    A > 0 || throw(ArgumentError("Cannot concatenate zero batches"))
+    return N, A
+end
+
+function DI.stack_batches_col(
+        block::F, ::DI.BatchSizeSettings{B, singlebatch, aligned}, ::NTuple{B, T}
+    ) where {F, B, singlebatch, aligned, T <: StaticArray}
+    N, A = static_batch_count(T, B)
+    M = hcat(ntuple(block, Val(A))...)
+    if aligned
+        return M
+    else
+        return M[:, SOneTo(N)]
+    end
+end
+
+function DI.stack_batches_row(
+        block::F, ::DI.BatchSizeSettings{B, singlebatch, aligned}, ::NTuple{B, T}
+    ) where {F, B, singlebatch, aligned, T <: StaticArray}
+    N, A = static_batch_count(T, B)
+    M = vcat(ntuple(block, Val(A))...)
+    if aligned
+        return M
+    else
+        return M[SOneTo(N), :]
+    end
 end
 
 DI.ismutable_array(::Type{<:SArray}) = false
diff --git a/DifferentiationInterface/src/first_order/jacobian.jl b/DifferentiationInterface/src/first_order/jacobian.jl
index 4c64b45bc..9e1cfc88e 100644
--- a/DifferentiationInterface/src/first_order/jacobian.jl
+++ b/DifferentiationInterface/src/first_order/jacobian.jl
@@ -354,22 +354,16 @@ end
 
 function _jacobian_aux(
         f_or_f!y::FY,
-        prep::PushforwardJacobianPrep{SIG, <:BatchSizeSettings{B, true, aligned}},
+        prep::PushforwardJacobianPrep{SIG, <:BatchSizeSettings{B, true}},
         backend::AbstractADType,
         x,
         contexts::Vararg{Context, C},
-    ) where {FY, SIG, B, aligned, C}
-    (; batch_size_settings, batched_seeds, pushforward_prep) = prep
-    (; B_last) = batch_size_settings
+    ) where {FY, SIG, B, C}
+    (; batched_seeds, pushforward_prep) = prep
     dy_batch = pushforward(
         f_or_f!y..., pushforward_prep, backend, x, only(batched_seeds), contexts...
     )
-    block = stack_vec_col(dy_batch)
-    if aligned
-        return block
-    else
-        return block[:, 1:B_last]
-    end
+    return stack_vec_col(dy_batch)
 end
 
 function _jacobian_aux(
@@ -379,13 +373,13 @@ function _jacobian_aux(
         x,
         contexts::Vararg{Context, C},
     ) where {FY, SIG, C}
-    (; batched_seeds, seed_example, pushforward_prep) = prep
+    (; batch_size_settings, batched_seeds, seed_example, pushforward_prep) = prep
 
     pushforward_prep_same = prepare_pushforward_same_point(
         f_or_f!y..., pushforward_prep, backend, x, seed_example, contexts...
     )
 
-    jac = stack(eachindex(batched_seeds); dims = 2) do a
+    return stack_batches_col(batch_size_settings, seed_example) do a
         dy = only(
             pushforward(
                 f_or_f!y...,
@@ -396,26 +390,24 @@ function _jacobian_aux(
                 contexts...,
             )
         )
-        return vec(dy)
+        return vec_block(dy)
     end
-    return jac
 end
 
 function _jacobian_aux(
         f_or_f!y::FY,
-        prep::PushforwardJacobianPrep{SIG, <:BatchSizeSettings{B, false, aligned}},
+        prep::PushforwardJacobianPrep{SIG, <:BatchSizeSettings{B, false}},
         backend::AbstractADType,
         x,
         contexts::Vararg{Context, C},
-    ) where {FY, SIG, B, aligned, C}
+    ) where {FY, SIG, B, C}
     (; batch_size_settings, batched_seeds, seed_example, pushforward_prep) = prep
-    (; A, B_last) = batch_size_settings
 
     pushforward_prep_same = prepare_pushforward_same_point(
         f_or_f!y..., pushforward_prep, backend, x, seed_example, contexts...
     )
 
-    jac = mapreduce(hcat, eachindex(batched_seeds)) do a
+    return stack_batches_col(batch_size_settings, seed_example) do a
         dy_batch = pushforward(
             f_or_f!y...,
             pushforward_prep_same,
@@ -424,37 +416,25 @@ function _jacobian_aux(
             batched_seeds[a],
             contexts...,
         )
-        block = stack_vec_col(dy_batch)
-        if !aligned && a == A
-            return block[:, 1:B_last]
-        else
-            return block
-        end
+        return stack_vec_col(dy_batch)
     end
-    return jac
 end
 
 function _jacobian_aux(
         f_or_f!y::FY,
-        prep::PullbackJacobianPrep{SIG, <:BatchSizeSettings{B, true, aligned}},
+        prep::PullbackJacobianPrep{SIG, <:BatchSizeSettings{B, true}},
         backend::AbstractADType,
         x,
         contexts::Vararg{Context, C},
-    ) where {FY, SIG, B, aligned, C}
-    (; batch_size_settings, batched_seeds, pullback_prep) = prep
-    (; B_last) = batch_size_settings
+    ) where {FY, SIG, B, C}
+    (; batched_seeds, pullback_prep) = prep
     dx_batch = pullback(
         f_or_f!y..., pullback_prep, backend, x, only(batched_seeds), contexts...
     )
     if eltype(x) <: Complex
         dx_batch = map(conj, dx_batch)
     end
-    block = stack_vec_row(dx_batch)
-    if aligned
-        return block
-    else
-        return block[1:B_last, :]
-    end
+    return stack_vec_row(dx_batch)
 end
 
 function _jacobian_aux(
@@ -464,56 +444,49 @@ function _jacobian_aux(
         x,
         contexts::Vararg{Context, C},
     ) where {FY, SIG, C}
-    (; batched_seeds, seed_example, pullback_prep) = prep
+    (; batch_size_settings, batched_seeds, seed_example, pullback_prep) = prep
 
     pullback_prep_same = prepare_pullback_same_point(
         f_or_f!y..., pullback_prep, backend, x, seed_example, contexts...
     )
 
-    jac = stack(eachindex(batched_seeds); dims = 1) do a
+    return stack_batches_row(batch_size_settings, seed_example) do a
         dx = only(
             pullback(
                 f_or_f!y..., pullback_prep_same, backend, x, batched_seeds[a], contexts...
             )
         )
+        # each block is a row of the Jacobian
         if eltype(x) <: Complex
-            return map(conj, vec(dx))
+            return transpose(map(conj, vec_block(dx)))
         else
-            return vec(dx)
+            return transpose(vec_block(dx))
         end
     end
-    return jac
 end
 
 function _jacobian_aux(
         f_or_f!y::FY,
-        prep::PullbackJacobianPrep{SIG, <:BatchSizeSettings{B, false, aligned}},
+        prep::PullbackJacobianPrep{SIG, <:BatchSizeSettings{B, false}},
         backend::AbstractADType,
         x,
         contexts::Vararg{Context, C},
-    ) where {FY, SIG, B, aligned, C}
+    ) where {FY, SIG, B, C}
     (; batch_size_settings, batched_seeds, seed_example, pullback_prep) = prep
-    (; A, B_last) = batch_size_settings
 
     pullback_prep_same = prepare_pullback_same_point(
         f_or_f!y..., pullback_prep, backend, x, seed_example, contexts...
     )
 
-    jac = mapreduce(vcat, eachindex(batched_seeds)) do a
+    return stack_batches_row(batch_size_settings, seed_example) do a
         dx_batch = pullback(
             f_or_f!y..., pullback_prep_same, backend, x, batched_seeds[a], contexts...
         )
         if eltype(x) <: Complex
             dx_batch = map(conj, dx_batch)
         end
-        block = stack_vec_row(dx_batch)
-        if !aligned && a == A
-            return block[1:B_last, :]
-        else
-            return block
-        end
+        return stack_vec_row(dx_batch)
     end
-    return jac
 end
 
 function _jacobian_aux!(
diff --git a/DifferentiationInterface/src/second_order/hessian.jl b/DifferentiationInterface/src/second_order/hessian.jl
index 84bbcad20..d2577acf5 100644
--- a/DifferentiationInterface/src/second_order/hessian.jl
+++ b/DifferentiationInterface/src/second_order/hessian.jl
@@ -159,44 +159,36 @@ function hessian(
         contexts::Vararg{Context, C},
     ) where {F, SIG, C}
     check_prep(f, prep, backend, x, contexts...)
-    (; batched_seeds, seed_example, hvp_prep) = prep
+    (; batch_size_settings, batched_seeds, seed_example, hvp_prep) = prep
 
     hvp_prep_same = prepare_hvp_same_point(
         f, hvp_prep, backend, x, seed_example, contexts...
     )
 
-    hess = mapreduce(hcat, eachindex(batched_seeds)) do a
+    return stack_batches_col(batch_size_settings, seed_example) do a
         dg = only(hvp(f, hvp_prep_same, backend, x, batched_seeds[a], contexts...))
-        return vec(dg)
+        return vec_block(dg)
     end
-    return hess
 end
 
 function hessian(
         f::F,
-        prep::HVPGradientHessianPrep{SIG, <:BatchSizeSettings{B, false, aligned}},
+        prep::HVPGradientHessianPrep{SIG, <:BatchSizeSettings{B, false}},
         backend::AbstractADType,
         x,
         contexts::Vararg{Context, C},
-    ) where {F, SIG, B, aligned, C}
+    ) where {F, SIG, B, C}
     check_prep(f, prep, backend, x, contexts...)
     (; batch_size_settings, batched_seeds, seed_example, hvp_prep) = prep
-    (; A, B_last) = batch_size_settings
 
     hvp_prep_same = prepare_hvp_same_point(
         f, hvp_prep, backend, x, seed_example, contexts...
     )
 
-    hess = mapreduce(hcat, eachindex(batched_seeds)) do a
+    return stack_batches_col(batch_size_settings, seed_example) do a
         dg_batch = hvp(f, hvp_prep_same, backend, x, batched_seeds[a], contexts...)
-        block = stack_vec_col(dg_batch)
-        if !aligned && a == A
-            return block[:, 1:B_last]
-        else
-            return block
-        end
+        return stack_vec_col(dg_batch)
     end
-    return hess
 end
 
 function hessian!(
diff --git a/DifferentiationInterface/src/utils/batchsize.jl b/DifferentiationInterface/src/utils/batchsize.jl
index 64c4288b4..20f8a03c6 100644
--- a/DifferentiationInterface/src/utils/batchsize.jl
+++ b/DifferentiationInterface/src/utils/batchsize.jl
@@ -13,23 +13,16 @@ Configuration for the batch size deduced from a backend and a sample array of le
 
   - `N::Int`: array length
   - `A::Int`: number of batches `A = div(N, B, RoundUp)`
-  - `B_last::Int`: size of the last batch (if `aligned` is `false`)
 """
 struct BatchSizeSettings{B, singlebatch, aligned}
     N::Int
     A::Int
-    B_last::Int
 end
 
 function BatchSizeSettings{B, singlebatch, aligned}(N::Integer) where {B, singlebatch, aligned}
     B > N > 0 && throw(ArgumentError("Batch size $B larger than input size $N"))
-    if B == N == 0
-        A = B_last = 0
-    else
-        A = div(N, B, RoundUp)
-        B_last = N % B
-    end
-    return BatchSizeSettings{B, singlebatch, aligned}(N, A, B_last)
+    A = B == N == 0 ? 0 : div(N, B, RoundUp)
+    return BatchSizeSettings{B, singlebatch, aligned}(N, A)
 end
 
 function BatchSizeSettings(B::Int, N::Int)
diff --git a/DifferentiationInterface/src/utils/linalg.jl b/DifferentiationInterface/src/utils/linalg.jl
index b27706a5e..85520fd71 100644
--- a/DifferentiationInterface/src/utils/linalg.jl
+++ b/DifferentiationInterface/src/utils/linalg.jl
@@ -1,6 +1,46 @@
+"""
+    vec_block(a)
+
+Flatten `a` into a vector suitable for concatenation into a Jacobian or Hessian block.
+"""
+vec_block(a) = vec(a)
+
 stack_vec_col(t::NTuple) = stack(vec, t; dims = 2)
 stack_vec_row(t::NTuple) = stack(vec, t; dims = 1)
 
+"""
+    stack_batches_col(block, batch_size_settings, seed_example)
+    stack_batches_row(block, batch_size_settings, seed_example)
+
+Concatenate the blocks `block(a)` computed from all batches `a` along columns (resp. rows), and drop the trailing entries of the last batch if it was padded by wrapping around the seeds.
+
+The `seed_example` is only used for dispatch, e.g. to deduce the number of seeds statically.
+
+Unlike `mapreduce(hcat, ...)`, `reduce(hcat, ...)` allocates the result only once.
+"""
+function stack_batches_col(
+        block::F, batch_size_settings::BatchSizeSettings{B, singlebatch, aligned}, _
+    ) where {F, B, singlebatch, aligned}
+    # a comprehension, since `map` fails to infer when `block` calls `map` itself
+    M = reduce(hcat, [block(a) for a in 1:(batch_size_settings.A)])
+    if aligned
+        return M
+    else
+        return M[:, 1:(batch_size_settings.N)]
+    end
+end
+
+function stack_batches_row(
+        block::F, batch_size_settings::BatchSizeSettings{B, singlebatch, aligned}, _
+    ) where {F, B, singlebatch, aligned}
+    M = reduce(vcat, [block(a) for a in 1:(batch_size_settings.A)])
+    if aligned
+        return M
+    else
+        return M[1:(batch_size_settings.N), :]
+    end
+end
+
 """
     ismutable_array(x)
 
diff --git a/DifferentiationInterface/test/Back/ForwardDiff/test.jl b/DifferentiationInterface/test/Back/ForwardDiff/test.jl
index c829e0555..e8e4b0bce 100644
--- a/DifferentiationInterface/test/Back/ForwardDiff/test.jl
+++ b/DifferentiationInterface/test/Back/ForwardDiff/test.jl
@@ -6,7 +6,7 @@ using DifferentiationInterface, DifferentiationInterfaceTest
 import DifferentiationInterface as DI
 import DifferentiationInterfaceTest as DIT
 using ForwardDiff: ForwardDiff
-using StaticArrays: StaticArrays, @SVector
+using StaticArrays: StaticArrays, @MVector, @SMatrix, @SVector, StaticMatrix
 using JLArrays: JLArrays
 using Test
 import JET, Chairmarks
@@ -94,6 +94,53 @@ end
             )
         ) isa DI.BatchSizeSettings{5}
     end
+
+    @testset "Static Jacobian with multiple batches" begin
+        f(x) = vcat(vec(x), vec(x) .^ 2)  # non-square Jacobian
+        function f!(y, x)
+            y[1:length(x)] .= vec(x)
+            y[(length(x) + 1):end] .= vec(x) .^ 2
+            return nothing
+        end
+        # a plain `AutoForwardDiff` calls `ForwardDiff.jacobian` and skips the batch assembly
+        @testset "$(typeof(x))" for x in (@SVector(rand(6)), @SMatrix(rand(2, 3)), @MVector(rand(6)))
+            y = StaticArrays.MVector{12, Float64}(undef)
+            J_ref = jacobian(f, AutoForwardDiff(), Array(x))
+            # B = 1, aligned, unaligned and single batch use different methods
+            @testset "Pushforward, chunksize $chunksize" for chunksize in (1, 3, 4, 6)
+                backend = DI.AutoForwardFromPrimitive(AutoForwardDiff(; chunksize))
+                J = @inferred jacobian(f, backend, x)
+                @test J isa StaticMatrix
+                @test J ≈ J_ref
+                J = @inferred jacobian(f!, y, backend, x)
+                @test J isa StaticMatrix
+                @test J ≈ J_ref
+            end
+            @testset "Pullback, chunksize $chunksize" for chunksize in (1, 4, 5, 12)
+                backend = DI.AutoReverseFromPrimitive(AutoForwardDiff(; chunksize))
+                J = @inferred jacobian(f, backend, x)
+                @test J isa StaticMatrix
+                @test J ≈ J_ref
+                J = @inferred jacobian(f!, y, backend, x)
+                @test J isa StaticMatrix
+                @test J ≈ J_ref
+            end
+        end
+    end
+
+    @testset "Static Hessian with multiple batches" begin
+        f(x) = sum(abs2, x) * sum(x)  # non-diagonal Hessian
+        for x in (@SVector(rand(6)), @SMatrix(rand(2, 3)), @MVector(rand(6)))
+            H_ref = hessian(f, AutoForwardDiff(), Array(x))
+            # B = 1, unaligned and single batch use different methods
+            @testset "Chunksize $chunksize" for chunksize in (1, 3, 4, 6)
+                backend = SecondOrder(DI.AutoSimpleFiniteDiff(; chunksize), AutoForwardDiff())
+                H = @inferred hessian(f, backend, x)
+                @test H isa StaticMatrix
+                @test H ≈ H_ref rtol = 1.0e-5
+            end
+        end
+    end
 end
 
 @testset verbose = true "Overloaded inputs" begin
diff --git a/DifferentiationInterface/test/Core/Internals/batchsize.jl b/DifferentiationInterface/test/Core/Internals/batchsize.jl
index 9a08aa5c1..49990e25c 100644
--- a/DifferentiationInterface/test/Core/Internals/batchsize.jl
+++ b/DifferentiationInterface/test/Core/Internals/batchsize.jl
@@ -7,6 +7,7 @@ using DifferentiationInterface:
     reasonable_batchsize,
     threshold_batchsize
 import DifferentiationInterface as DI
+using LinearAlgebra: I
 using StaticArrays
 using Test
 
@@ -26,12 +27,12 @@ BSS = BatchSizeSettings
 end
 
 @testset "Constructor" begin
-    @test BSS(1, 0) === BSS{1, false, true}(0, 0, 0)
-    @test BSS(0, 0) === BSS{0, true, true}(0, 0, 0)
-    @test BSS(1, 1) === BSS{1, true, true}(1, 1, 0)
-    @test BSS(2, 2) === BSS{2, true, true}(2, 1, 0)
-    @test BSS(4, 8) === BSS{4, false, true}(8, 2, 0)
-    @test BSS(8, 13) === BSS{8, false, false}(13, 2, 5)
+    @test BSS(1, 0) === BSS{1, false, true}(0, 0)
+    @test BSS(0, 0) === BSS{0, true, true}(0, 0)
+    @test BSS(1, 1) === BSS{1, true, true}(1, 1)
+    @test BSS(2, 2) === BSS{2, true, true}(2, 1)
+    @test BSS(4, 8) === BSS{4, false, true}(8, 2)
+    @test BSS(8, 13) === BSS{8, false, false}(13, 2)
     @test_throws ArgumentError BSS(3, 2)
 end
 
@@ -87,12 +88,71 @@ end
             @inferred prepare_jacobian(identity, AutoSimpleFiniteDiff(; chunksize = 4), x)
         ) isa DI.PushforwardJacobianPrep{<:Any, BSS{4, false, false}}
         @test (@inferred jacobian(identity, AutoSimpleFiniteDiff(), x)) isa StaticMatrix
-        # multiple batches are assembled with `mapreduce(hcat, ...)`, which changes the
-        # static size at every step
-        @test_broken (
-            @inferred jacobian(identity, AutoSimpleFiniteDiff(; chunksize = 4), x)
-        ) isa StaticMatrix
     end
+
+    # multiple batches must be concatenated without losing static sizes or inference
+    f(x) = vcat(vec(x), vec(x) .^ 2)  # non-square Jacobian
+    function f!(y, x)
+        y[1:length(x)] .= vec(x)
+        y[(length(x) + 1):end] .= vec(x) .^ 2
+        return nothing
+    end
+    g(x) = sum(abs2, x) * sum(x)  # non-diagonal Hessian
+    has_type(J, x) = x isa StaticArray ? J isa Union{SArray, MArray} : J isa Matrix
+    # without preparation, the batch size is only inferred for static arrays
+    function inferred_jacobian(f_or_f!y, backend, x)
+        x isa StaticArray && @inferred jacobian(f_or_f!y..., backend, x)
+        prep = prepare_jacobian(f_or_f!y..., backend, x)
+        return @inferred jacobian(f_or_f!y..., prep, backend, x)
+    end
+    function inferred_hessian(f, backend, x)
+        x isa StaticArray && @inferred hessian(f, backend, x)
+        prep = prepare_hessian(f, backend, x)
+        return @inferred hessian(f, prep, backend, x)
+    end
+    @testset "$(typeof(x))" for x in (
+            @SVector(rand(6)),
+            @SMatrix(rand(2, 3)),
+            @MVector(rand(6)),
+            rand(6),
+            @SVector(rand(ComplexF64, 6)),
+        )
+        y = x isa StaticArray ? MVector{12, eltype(x)}(undef) : similar(x, 12)
+        J_ref = jacobian(f, AutoSimpleFiniteDiff(), Array(x))
+        # B = 1, aligned, unaligned and single batch use different methods
+        @testset "Pushforward, chunksize $chunksize" for chunksize in (1, 3, 4, 6)
+            backend = AutoSimpleFiniteDiff(; chunksize)
+            J = inferred_jacobian((f,), backend, x)
+            @test has_type(J, x)
+            @test J ≈ J_ref
+            J = inferred_jacobian((f!, y), backend, x)
+            @test has_type(J, x)
+            @test J ≈ J_ref
+        end
+        @testset "Pullback, chunksize $chunksize" for chunksize in (1, 4, 5, 12)
+            backend = DI.AutoReverseFromPrimitive(AutoSimpleFiniteDiff(; chunksize))
+            J = inferred_jacobian((f,), backend, x)
+            @test has_type(J, x)
+            @test J ≈ J_ref
+            J = inferred_jacobian((f!, y), backend, x)
+            @test has_type(J, x)
+            @test J ≈ J_ref
+        end
+        if eltype(x) <: Real
+            H_ref = 2 * (sum(x) * I + (vec(x) .+ vec(x)'))
+            @testset "Hessian, chunksize $chunksize" for chunksize in (1, 3, 4, 6)
+                backend = SecondOrder(
+                    AutoSimpleFiniteDiff(; chunksize), AutoSimpleFiniteDiff()
+                )
+                H = inferred_hessian(g, backend, x)
+                @test has_type(H, x)
+                @test H ≈ H_ref rtol = 1.0e-4
+            end
+        end
+    end
+    @test_throws ArgumentError("Cannot concatenate zero batches") jacobian(
+        identity, AutoSimpleFiniteDiff(; chunksize = 2), @SVector(zeros(0))
+    )
 end
 
 @testset "Thresholding" begin
diff --git a/DifferentiationInterface/test/Core/SimpleFiniteDiff/test.jl b/DifferentiationInterface/test/Core/SimpleFiniteDiff/test.jl
index 1ff483e17..05e9da732 100644
--- a/DifferentiationInterface/test/Core/SimpleFiniteDiff/test.jl
+++ b/DifferentiationInterface/test/Core/SimpleFiniteDiff/test.jl
@@ -1,4 +1,5 @@
 using DifferentiationInterface, DifferentiationInterfaceTest
+import DifferentiationInterfaceTest as DIT
 using DifferentiationInterface:
     AutoSimpleFiniteDiff,
     AutoForwardFromPrimitive,
@@ -6,6 +7,7 @@ using DifferentiationInterface:
     DenseSparsityDetector
 using SparseMatrixColorings
 using JLArrays, StaticArrays
+using JET: JET
 using Test
 
 backends = [ #
@@ -168,6 +170,50 @@ end
         vcat(static_scenarios(), gpu_scenarios());
         logging = LOGGING,
     )
+    # multiple batches
+    forward_backend = AutoSimpleFiniteDiff(; chunksize = 2)
+    reverse_backend = AutoReverseFromPrimitive(AutoSimpleFiniteDiff(; chunksize = 2))
+    multibatch_scenarios = filter(
+        s -> DIT.operator(s) in (:jacobian, :hessian), static_scenarios()
+    )
+    test_differentiation(
+        [forward_backend, reverse_backend], multibatch_scenarios; logging = LOGGING
+    )
+
+    # JET < 0.12 does not define `JET_AVAILABLE`
+    if !isdefined(JET, :JET_AVAILABLE) || JET.JET_AVAILABLE
+        inplace_hessian(s) = DIT.operator(s) == :hessian && DIT.operator_place(s) == :in
+        test_differentiation(
+            forward_backend,
+            multibatch_scenarios;
+            correctness = false,
+            type_stability = :prepared,
+            logging = LOGGING,
+        )
+        test_differentiation(
+            reverse_backend,
+            filter(!inplace_hessian, multibatch_scenarios);
+            correctness = false,
+            type_stability = :prepared,
+            logging = LOGGING,
+        )
+        @testset "Type stability (broken)" for scen in filter(inplace_hessian, multibatch_scenarios)
+            (; f, x, res1, res2, contexts, prep_args) = deepcopy(scen)
+            prep = prepare_hessian(f, reverse_backend, prep_args.x, prep_args.contexts...)
+            JET.@test_opt broken = true hessian!(
+                f, DIT.mysimilar(res2), prep, reverse_backend, x, contexts...
+            )
+            JET.@test_opt broken = true value_gradient_and_hessian!(
+                f,
+                DIT.mysimilar(res1),
+                DIT.mysimilar(res2),
+                prep,
+                reverse_backend,
+                x,
+                contexts...,
+            )
+        end
+    end
 end;
 
 @testset "Array format preservation in wrong mode" begin