Allocations for non-adapating BVP.solve() using MIRK

[maxschelski]

I posted this first as a question in a Julia Discourse thread:
https://discourse.julialang.org/t/boundaryvaluediffeq-jl-reducing-allocations/136255

I noticed quite a few allocations when using a MIRK solver and turning adaptivity off.
I turn adaptivity off for my BVP, since I seem to have to set the grid size manually to get accurate results.

Minimum reproducible example

import BoundaryValueDiffEq as BVP


function distribution_fun_test!(dydt::AbstractVector{T}, y::AbstractVector{T},  p::Any, #
                                t::Float64) where {T<:Real}

    y1, y2, y3, y4 = y

    dydt[1] = 0.001 * y2
    dydt[2] = 10 * y1
    dydt[3] = -0.03 * y2 + 0.0001 * y4
    dydt[4] = 10 * y3
                    
    uptake = 0.1 * y3 / (1000 * y3 + 1)
    dydt[3] += uptake

end


function distribution_bc_test!(residual::AbstractVector{}, u::Any, p::Any, t::Vector{Float64})

    residual[1] = (u[4, end] - 0.001)
    residual[2] = u[1, end]
    residual[3] =   (- u[3, end])
    residual[4] =   (u[1, 1])

end

function test_solve()
    L = 100.0
    mesh_size = 500
    tspan::Tuple{Float64, Float64} = (0.0, Float64(L))

    y_initial::Vector{Float64} = [1, 0.001/0.03, 1, 0.001]

    bvp_problem = BVP.BVProblem(distribution_fun_test!, 
                                distribution_bc_test!, 
                                y_initial, tspan)
        
    @time BVP.solve(bvp_problem, BVP.MIRK5(); adaptive = false, dt=Float64(L/mesh_size))
    @time BVP.solve(bvp_problem, BVP.MIRK5(); adaptive = false, dt=Float64(L/mesh_size))
end

test_solve()

Observation
The output of the timing I am getting is:

3.606669 seconds (4.35 M allocations: 216.794 MiB, 1.00% gc time, 96.31% compilation time: 100% of which was recompilation) 
0.132779 seconds (440.05 k allocations: 26.994 MiB)

Which shows 440k allocations in the second iteration, with adaptivity turned off.
As Chris Rackauckas suggested in the thread above, maybe it could be tested whether making the solve without adaptivity non-allocating would improve performance?

Environment
Windows 10
Julia v1.12.
BoundaryValueDiffEq v5.18.0