diff --git a/benchmark/simple_pendulum.jl b/benchmark/simple_pendulum.jl
index e2151a71..8663fbed 100644
--- a/benchmark/simple_pendulum.jl
+++ b/benchmark/simple_pendulum.jl
@@ -18,12 +18,6 @@ function bc_pendulum!(residual, u, p, t)
     return nothing
 end
 
-function bc_pendulum_mirk!(residual, u, p, t)
-    residual[1] = u[:, end ÷ 2][1] + π / 2
-    residual[2] = u[:, end][1] - π / 2
-    return nothing
-end
-
 function simple_pendulum(u, p, t)
     g, L, θ, dθ = 9.81, 1.0, u[1], u[2]
     return [dθ, -(g / L) * sin(θ)]
@@ -34,16 +28,8 @@ function bc_pendulum(u, p, t)
     return [u((t0 + t1) / 2)[1] + π / 2, u(t1)[1] - π / 2]
 end
 
-function bc_pendulum_mirk(u, p, t)
-    return [u[:, end ÷ 2][1] + π / 2, u[:, end][1] - π / 2]
-end
-
 const prob_oop = BVProblem{false}(simple_pendulum, bc_pendulum, [π / 2, π / 2], tspan)
 const prob_iip = BVProblem{true}(simple_pendulum!, bc_pendulum!, [π / 2, π / 2], tspan)
-const prob_oop_mirk = BVProblem{false}(
-    simple_pendulum, bc_pendulum_mirk, [π / 2, π / 2], tspan)
-const prob_iip_mirk = BVProblem{true}(
-    simple_pendulum!, bc_pendulum_mirk!, [π / 2, π / 2], tspan)
 
 end
 
@@ -77,7 +63,7 @@ function create_simple_pendulum_benchmark()
     for alg in (MIRK2, MIRK3, MIRK4, MIRK5, MIRK6)
         if @isdefined(alg)
             iip_suite["$alg()"] = @benchmarkable solve(
-                $SimplePendulumBenchmark.prob_iip_mirk, $alg(), dt = 0.05)
+                $SimplePendulumBenchmark.prob_iip, $alg(), dt = 0.05)
         end
     end
 
@@ -102,7 +88,7 @@ function create_simple_pendulum_benchmark()
     for alg in (MIRK2, MIRK3, MIRK4, MIRK5, MIRK6)
         if @isdefined(alg)
             oop_suite["$alg()"] = @benchmarkable solve(
-                $SimplePendulumBenchmark.prob_oop_mirk, $alg(), dt = 0.05)
+                $SimplePendulumBenchmark.prob_oop, $alg(), dt = 0.05)
         end
     end
 
diff --git a/lib/BoundaryValueDiffEqCore/src/misc_utils.jl b/lib/BoundaryValueDiffEqCore/src/misc_utils.jl
index fed90238..4f06e27e 100644
--- a/lib/BoundaryValueDiffEqCore/src/misc_utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/misc_utils.jl
@@ -1,7 +1,14 @@
 # Intermidiate solution evaluation
-@concrete struct EvalSol{iip}
+@concrete struct EvalSol{C}
     u
     t
-    alg
-    k_discrete
+    cache::C
 end
+
+Base.size(e::EvalSol) = (size(e.u[1])..., length(e.u))
+Base.size(e::EvalSol, i) = size(e)[i]
+
+Base.axes(e::EvalSol) = Base.OneTo.(size(e))
+Base.axes(e::EvalSol, d::Int) = Base.OneTo.(size(e)[d])
+
+Base.getindex(e::EvalSol, args...) = Base.getindex(VectorOfArray(e.u), args...)
diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 41d07abc..33afd602 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -223,7 +223,7 @@ end
 __vec_f(u, p, t, f, u_size) = vec(f(reshape(u, u_size), p, t))
 
 function __vec_bc!(resid, sol, p, t, bc!, resid_size, u_size)
-    bc!(reshape(resid, resid_size), __restructure_sol(sol, u_size), p, t)
+    bc!(reshape(resid, resid_size), sol, p, t)
     return nothing
 end
 
@@ -232,7 +232,7 @@ function __vec_bc!(resid, sol, p, bc!, resid_size, u_size)
     return nothing
 end
 
-__vec_bc(sol, p, t, bc, u_size) = vec(bc(__restructure_sol(sol, u_size), p, t))
+__vec_bc(sol, p, t, bc, u_size) = vec(bc(sol, p, t))
 __vec_bc(sol, p, bc, u_size) = vec(bc(reshape(sol, u_size), p))
 
 @inline __get_non_sparse_ad(ad::AbstractADType) = ad
@@ -240,9 +240,11 @@ __vec_bc(sol, p, bc, u_size) = vec(bc(reshape(sol, u_size), p))
 
 # Restructure Solution
 function __restructure_sol(sol::AbstractVectorOfArray, u_size)
+    (size(first(sol)) == u_size) && return sol
     return VectorOfArray(map(Base.Fix2(reshape, u_size), sol))
 end
-function __restructure_sol(sol::Vector{<:AbstractArray}, u_size)
+function __restructure_sol(sol::AbstractArray{<:AbstractArray}, u_size)
+    (size(first(sol)) == u_size) && return sol
     return map(Base.Fix2(reshape, u_size), sol)
 end
 
diff --git a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
index 4cb76c9b..860a3af9 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
@@ -15,15 +15,16 @@ import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorit
                                 recursive_flatten, recursive_flatten!, recursive_unflatten!,
                                 __concrete_nonlinearsolve_algorithm, diff!,
                                 __FastShortcutBVPCompatibleNonlinearPolyalg,
-                                __FastShortcutBVPCompatibleNLLSPolyalg,
+                                __FastShortcutBVPCompatibleNLLSPolyalg, EvalSol,
                                 concrete_jacobian_algorithm, eval_bc_residual,
                                 eval_bc_residual!, get_tmp, __maybe_matmul!,
                                 __append_similar!, __extract_problem_details,
                                 __initial_guess, __maybe_allocate_diffcache,
-                                __get_bcresid_prototype, __similar, __vec, __vec_f,
-                                __vec_f!, __vec_bc, __vec_bc!, recursive_flatten_twopoint!,
-                                __internal_nlsolve_problem, MaybeDiffCache, __extract_mesh,
-                                __extract_u0, __has_initial_guess, __initial_guess_length,
+                                __restructure_sol, __get_bcresid_prototype, __similar,
+                                __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
+                                recursive_flatten_twopoint!, __internal_nlsolve_problem,
+                                MaybeDiffCache, __extract_mesh, __extract_u0,
+                                __has_initial_guess, __initial_guess_length,
                                 __initial_guess_on_mesh, __flatten_initial_guess,
                                 __build_solution, __Fix3, __sparse_jacobian_cache,
                                 __sparsity_detection_alg, _sparse_like, ColoredMatrix
@@ -33,7 +34,7 @@ import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scal
 import ConcreteStructs: @concrete
 import DiffEqBase: solve
 import FastClosures: @closure
-import ForwardDiff: ForwardDiff, pickchunksize
+import ForwardDiff: ForwardDiff, pickchunksize, Dual
 import Logging
 import RecursiveArrayTools: ArrayPartition, DiffEqArray
 import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
@@ -60,11 +61,11 @@ include("sparse_jacobians.jl")
     f1 = (u, p, t) -> [u[2], 0]
 
     function bc1!(residual, u, p, t)
-        residual[1] = u[:, 1][1] - 5
-        residual[2] = u[:, end][1]
+        residual[1] = u(0.0)[1] - 5
+        residual[2] = u(5.0)[1]
     end
 
-    bc1 = (u, p, t) -> [u[:, 1][1] - 5, u[:, end][1]]
+    bc1 = (u, p, t) -> [u(0.0)[1] - 5, u(5.0)[1]]
 
     bc1_a! = (residual, ua, p) -> (residual[1] = ua[1] - 5)
     bc1_b! = (residual, ub, p) -> (residual[1] = ub[1])
@@ -143,14 +144,14 @@ include("sparse_jacobians.jl")
     f1_nlls = (u, p, t) -> [u[2], -u[1]]
 
     bc1_nlls! = (resid, sol, p, t) -> begin
-        solₜ₁ = sol[:, 1]
-        solₜ₂ = sol[:, end]
+        solₜ₁ = sol(0.0)
+        solₜ₂ = sol(100.0)
         resid[1] = solₜ₁[1]
         resid[2] = solₜ₂[1] - 1
         resid[3] = solₜ₂[2] + 1.729109
         return nothing
     end
-    bc1_nlls = (sol, p, t) -> [sol[:, 1][1], sol[:, end][1] - 1, sol[:, end][2] + 1.729109]
+    bc1_nlls = (sol, p, t) -> [sol(0.0)[1], sol(100.0)[1] - 1, sol(100.0)[2] + 1.729109]
 
     bc1_nlls_a! = (resid, ua, p) -> (resid[1] = ua[1])
     bc1_nlls_b! = (resid, ub, p) -> (resid[1] = ub[1] - 1;
diff --git a/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl b/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl
index 033e942e..014101c8 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl
@@ -390,10 +390,10 @@ function apply_q_prime(τ, h, coeffs)
     return sum(i * coeffs[i] * (τ * h)^(i - 1) for i in axes(coeffs, 1))
 end
 
-function eval_q(y_i, τ, h, A, K)
+function eval_q(y_i::AbstractArray{T}, τ, h, A, K) where {T}
     M = size(K, 1)
-    q = zeros(M)
-    q′ = zeros(M)
+    q = zeros(T, M)
+    q′ = zeros(T, M)
     for i in 1:M
         ki = @view K[i, :]
         coeffs = get_q_coeffs(A, ki, h)
diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 4cc2649b..404a9dae 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -353,7 +353,7 @@ end
 
 function __perform_firk_iteration(cache::Union{FIRKCacheExpand, FIRKCacheNested}, abstol,
         adaptive::Bool; nlsolve_kwargs = (;), kwargs...)
-    nlprob = __construct_nlproblem(cache, vec(cache.y₀))
+    nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
     nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve)
     sol_nlprob = __solve(
         nlprob, nlsolve_alg; abstol, kwargs..., nlsolve_kwargs..., alias_u0 = true)
@@ -402,9 +402,11 @@ end
 
 # Constructing the Nonlinear Problem
 function __construct_nlproblem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpand{iip}},
-        y::AbstractVector) where {iip}
+        y::AbstractVector, y₀::AbstractVectorOfArray) where {iip}
     pt = cache.problem_type
 
+    eval_sol = EvalSol(__restructure_sol(y₀.u, cache.in_size), cache.mesh, cache)
+
     loss_bc = if iip
         @closure (du, u, p) -> __firk_loss_bc!(
             du, u, p, pt, cache.bc, cache.y, cache.mesh, cache)
@@ -422,9 +424,10 @@ function __construct_nlproblem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpan
 
     loss = if iip
         @closure (du, u, p) -> __firk_loss!(
-            du, u, p, cache.y, pt, cache.bc, cache.residual, cache.mesh, cache)
+            du, u, p, cache.y, pt, cache.bc, cache.residual, cache.mesh, cache, eval_sol)
     else
-        @closure (u, p) -> __firk_loss(u, p, cache.y, pt, cache.bc, cache.mesh, cache)
+        @closure (u, p) -> __firk_loss(
+            u, p, cache.y, pt, cache.bc, cache.mesh, cache, eval_sol)
     end
 
     return __construct_nlproblem(cache, y, loss_bc, loss_collocation, loss, pt)
@@ -658,19 +661,19 @@ function __construct_nlproblem(
     return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
 end
 
-@views function __firk_loss!(
-        resid, u, p, y, pt::StandardBVProblem, bc!::BC, residual, mesh, cache) where {BC}
+@views function __firk_loss!(resid, u, p, y, pt::StandardBVProblem, bc!::BC,
+        residual, mesh, cache, eval_sol) where {BC}
     y_ = recursive_unflatten!(y, u)
     resids = [get_tmp(r, u) for r in residual]
-    soly_ = VectorOfArray(y_)
-    eval_bc_residual!(resids[1], pt, bc!, soly_, p, mesh)
     Φ!(resids[2:end], cache, y_, u, p)
+    eval_sol.u[1:end] .= y_
+    eval_bc_residual!(resids[1], pt, bc!, eval_sol, p, mesh)
     recursive_flatten!(resid, resids)
     return nothing
 end
 
 @views function __firk_loss!(resid, u, p, y, pt::TwoPointBVProblem, bc!::Tuple{BC1, BC2},
-        residual, mesh, cache) where {BC1, BC2}
+        residual, mesh, cache, _) where {BC1, BC2}
     y_ = recursive_unflatten!(y, u)
     soly_ = VectorOfArray(y_)
     resids = [get_tmp(r, u) for r in residual]
@@ -682,16 +685,17 @@ end
     return nothing
 end
 
-@views function __firk_loss(u, p, y, pt::StandardBVProblem, bc::BC, mesh, cache) where {BC}
+@views function __firk_loss(
+        u, p, y, pt::StandardBVProblem, bc::BC, mesh, cache, eval_sol) where {BC}
     y_ = recursive_unflatten!(y, u)
-    soly_ = VectorOfArray(y_)
-    resid_bc = eval_bc_residual(pt, bc, soly_, p, mesh)
+    eval_sol.u[1:end] .= y_
+    resid_bc = eval_bc_residual(pt, bc, eval_sol, p, mesh)
     resid_co = Φ(cache, y_, u, p)
     return vcat(resid_bc, mapreduce(vec, vcat, resid_co))
 end
 
-@views function __firk_loss(
-        u, p, y, pt::TwoPointBVProblem, bc::Tuple{BC1, BC2}, mesh, cache) where {BC1, BC2}
+@views function __firk_loss(u, p, y, pt::TwoPointBVProblem, bc::Tuple{BC1, BC2},
+        mesh, cache, _) where {BC1, BC2}
     y_ = recursive_unflatten!(y, u)
     soly_ = VectorOfArray(y_)
     resid_bca, resid_bcb = eval_bc_residual(pt, bc, soly_, p, mesh)
@@ -702,16 +706,16 @@ end
 @views function __firk_loss_bc!(resid, u, p, pt, bc!::BC, y, mesh,
         cache::Union{FIRKCacheNested, FIRKCacheExpand}) where {BC}
     y_ = recursive_unflatten!(y, u)
-    soly_ = VectorOfArray(y_)
-    eval_bc_residual!(resid, pt, bc!, soly_, p, mesh)
+    eval_sol = EvalSol(__restructure_sol(y_, cache.in_size), mesh, cache)
+    eval_bc_residual!(resid, pt, bc!, eval_sol, p, mesh)
     return nothing
 end
 
 @views function __firk_loss_bc(u, p, pt, bc!::BC, y, mesh,
         cache::Union{FIRKCacheNested, FIRKCacheExpand}) where {BC}
     y_ = recursive_unflatten!(y, u)
-    soly_ = VectorOfArray(y_)
-    return eval_bc_residual(pt, bc!, soly_, p, mesh)
+    eval_sol = EvalSol(__restructure_sol(y_, cache.in_size), mesh, cache)
+    return eval_bc_residual(pt, bc!, eval_sol, p, mesh)
 end
 
 @views function __firk_loss_collocation!(resid, u, p, y, mesh, residual, cache)
diff --git a/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl
index fac00915..82eae516 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl
@@ -127,3 +127,97 @@ end
     cache.mesh, u, cache)
 @inline __build_interpolation(cache::FIRKCacheNested, u::AbstractVector) = FIRKNestedInterpolation(
     cache.mesh, u, cache)
+
+# Intermidiate solution for evaluating boundry conditions
+# basically simplified version of the interpolation for FIRK
+# Expanded FIRK
+function (s::EvalSol{C})(tval::Number) where {C <: FIRKCacheExpand}
+    (; t, u, cache) = s
+    (; f, alg, ITU, p) = cache
+    (; q_coeff) = ITU
+    stage = alg_stage(alg)
+    # Quick handle for the case where tval is at the boundary
+    (tval == t[1]) && return first(u)
+    (tval == t[end]) && return last(u)
+    K = __similar(first(u), length(first(u)), stage)
+    j = interval(t, tval)
+    ctr_y = (j - 1) * (stage + 1) + 1
+
+    yᵢ = u[ctr_y]
+    yᵢ₊₁ = u[ctr_y + stage + 1]
+
+    if SciMLBase.isinplace(cache.prob)
+        dyᵢ = similar(yᵢ)
+        dyᵢ₊₁ = similar(yᵢ₊₁)
+
+        f(dyᵢ, yᵢ, p, t[j])
+        f(dyᵢ₊₁, yᵢ₊₁, p, t[j + 1])
+    else
+        dyᵢ = f(yᵢ, p, t[j])
+        dyᵢ₊₁ = f(yᵢ₊₁, p, t[j + 1])
+    end
+
+    # Load interpolation residual
+    for jj in 1:stage
+        K[:, jj] = u[ctr_y + jj]
+    end
+    h = t[j + 1] - t[j]
+    τ = tval - t[j]
+
+    z₁, z₁′ = eval_q(yᵢ, 0.5, h, q_coeff, K) # Evaluate q(x) at midpoints
+    S_coeffs = get_S_coeffs(h, yᵢ, yᵢ₊₁, z₁, dyᵢ, dyᵢ₊₁, z₁′)
+
+    z = similar(yᵢ)
+
+    S_interpolate!(z, τ, S_coeffs)
+    return z
+end
+
+nodual_value(x) = x
+nodual_value(x::Dual) = ForwardDiff.value(x)
+nodual_value(x::AbstractArray{<:Dual}) = map(ForwardDiff.value, x)
+
+# Nested FIRK
+function (s::EvalSol{C})(tval::Number) where {C <: FIRKCacheNested}
+    (; t, u, cache) = s
+    (; f, nest_prob, nest_tol, alg, mesh_dt, p, ITU) = cache
+    (; q_coeff) = ITU
+    stage = alg_stage(alg)
+    # Quick handle for the case where tval is at the boundary
+    (tval == t[1]) && return first(u)
+    (tval == t[end]) && return last(u)
+    j = interval(t, tval)
+    h = t[j + 1] - t[j]
+    τ = tval - t[j]
+
+    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nestprob_p = zeros(cache.M + 2)
+
+    yᵢ = u[j]
+    yᵢ₊₁ = u[j + 1]
+
+    if SciMLBase.isinplace(cache.prob)
+        dyᵢ = similar(yᵢ)
+        dyᵢ₊₁ = similar(yᵢ₊₁)
+
+        f(dyᵢ, yᵢ, p, t[j])
+        f(dyᵢ₊₁, yᵢ₊₁, p, t[j + 1])
+    else
+        dyᵢ = f(yᵢ, p, t[j])
+        dyᵢ₊₁ = f(yᵢ₊₁, p, t[j + 1])
+    end
+
+    nestprob_p[1] = t[j]
+    nestprob_p[2] = mesh_dt[j]
+    nestprob_p[3:end] .= nodual_value(yᵢ)
+
+    _nestprob = remake(nest_prob, p = nestprob_p)
+    nestsol = __solve(_nestprob, nest_nlsolve_alg; abstol = nest_tol)
+    K = nestsol.u
+
+    z₁, z₁′ = eval_q(yᵢ, 0.5, h, q_coeff, K) # Evaluate q(x) at midpoints
+    S_coeffs = get_S_coeffs(h, yᵢ, yᵢ₊₁, z₁, dyᵢ, dyᵢ₊₁, z₁′)
+    z = similar(yᵢ)
+    S_interpolate!(z, τ, S_coeffs)
+    return z
+end
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
index 4f106445..d9f5623f 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
@@ -7,8 +7,8 @@
     end
 
     function bc!(residual, u, p, t)
-        residual[1] = u[:, 1][1] - 1.0
-        residual[2] = u[:, end][1]
+        residual[1] = u(0.0)[1] - 1.0
+        residual[2] = u(1.0)[1]
     end
 
     prob_func(prob, i, repeat) = remake(prob, p = [rand()])
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index fb5e9213..5876991c 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -39,10 +39,17 @@ end
 f2(u, p, t) = [u[2], -u[1]]
 
 function boundary!(residual, u, p, t)
+    residual[1] = u(0.0)[1] - 5
+    residual[2] = u(5.0)[1]
+end
+boundary(u, p, t) = [u(0.0)[1] - 5, u(5.0)[1]]
+
+# Array indexing for boudnary conditions
+function boundary_indexing!(residual, u, p, t)
     residual[1] = u[:, 1][1] - 5
     residual[2] = u[:, end][1]
 end
-boundary(u, p, t) = [u[:, 1][1] - 5, u[:, end][1]]
+boundary_indexing(u, p, t) = [u[:, 1][1] - 5, u[:, end][1]]
 
 function boundary_two_point_a!(resida, ua, p)
     resida[1] = ua[1] - 5
@@ -75,6 +82,8 @@ probArr = [BVProblem(odef1!, boundary!, u0, tspan, nlls = Val(false)),
     BVProblem(odef1, boundary, u0, tspan, nlls = Val(false)),
     BVProblem(odef2!, boundary!, u0, tspan, nlls = Val(false)),
     BVProblem(odef2, boundary, u0, tspan, nlls = Val(false)),
+    BVProblem(odef2!, boundary_indexing!, u0, tspan, nlls = Val(false)),
+    BVProblem(odef2, boundary_indexing, u0, tspan, nlls = Val(false)),
     TwoPointBVProblem(odef1!, (boundary_two_point_a!, boundary_two_point_b!),
         u0, tspan; bcresid_prototype, nlls = Val(false)),
     TwoPointBVProblem(odef1, (boundary_two_point_a, boundary_two_point_b),
@@ -96,7 +105,7 @@ end
 @testitem "Affineness" setup=[FIRKExpandedConvergenceTests] begin
     using LinearAlgebra
 
-    @testset "Problem: $i" for i in (1, 2, 5, 6)
+    @testset "Problem: $i" for i in (1, 2, 7, 8)
         prob = probArr[i]
 
         @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5)
@@ -125,7 +134,7 @@ end
 @testitem "JET: Runtime Dispatches" setup=[FIRKExpandedConvergenceTests] begin
     using JET
 
-    @testset "Problem: $i" for i in 1:8
+    @testset "Problem: $i" for i in 1:10
         prob = probArr[i]
         @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5)
             solver = lobattoIIIa_solver(Val(stage); nlsolve = NewtonRaphson(),
@@ -165,13 +174,13 @@ end
 @testitem "Convergence on Linear" setup=[FIRKExpandedConvergenceTests] begin
     using LinearAlgebra, DiffEqDevTools
 
-    @testset "Problem: $i" for i in (3, 4, 7, 8)
+    @testset "Problem: $i" for i in (3, 4, 9, 10)
         prob = probArr[i]
 
         @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5)
             @time sim = test_convergence(
                 dts, prob, lobattoIIIa_solver(Val(stage)); abstol = 1e-8)
-            if (stage == 5) || (((i == 7) || (i == 8)) && stage == 4)
+            if (stage == 5) || (((i == 9) || (i == 10)) && stage == 4)
                 @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             else
                 @test sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
@@ -181,7 +190,7 @@ end
         @testset "LobattoIIIb$stage" for stage in (2, 3, 4, 5)
             @time sim = test_convergence(
                 dts, prob, lobattoIIIb_solver(Val(stage)); abstol = 1e-8, reltol = 1e-8)
-            if (stage == 5) || (stage == 4 && i == 8)
+            if (stage == 5) || (stage == 4 && i == 10)
                 @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             elseif stage == 4
                 @test sim.𝒪est[:final]≈2 * stage - 2 atol=0.5
@@ -226,8 +235,8 @@ end
     end
 
     function bc_pendulum!(residual, u, p, t)
-        residual[1] = u[:, end ÷ 2][1] + π / 2 # the solution at the middle of the time span should be -pi/2
-        residual[2] = u[:, end][1] - π / 2 # the solution at the end of the time span should be pi/2
+        residual[1] = u(pi / 4)[1] + π / 2 # the solution at the middle of the time span should be -pi/2
+        residual[2] = u(pi / 2)[1] - π / 2 # the solution at the end of the time span should be pi/2
     end
 
     u0 = MVector{2}([pi / 2, pi / 2])
@@ -279,8 +288,8 @@ end
         du[2] = 1 / p * u[1]
     end
     function prob_bvp_linear_bc!(res, u, p, t)
-        res[1] = u[:, 1][1] - 1
-        res[2] = u[:, end][1]
+        res[1] = u(0.0)[1] - 1
+        res[2] = u(1.0)[1]
     end
 
     prob_bvp_linear_function = ODEFunction(
@@ -358,12 +367,12 @@ end
     end
 
     function swirling_flow_bc!(res, u, p, t)
-        res[1] = u[:, 1][1] + 1.0
-        res[2] = u[:, 1][3]
-        res[3] = u[:, 1][4]
-        res[4] = u[:, end][1] - 1.0
-        res[5] = u[:, end][3]
-        res[6] = u[:, end][4]
+        res[1] = u(0.0)[1] + 1.0
+        res[2] = u(0.0)[3]
+        res[3] = u(0.0)[4]
+        res[4] = u(1.0)[1] - 1.0
+        res[5] = u(1.0)[3]
+        res[6] = u(1.0)[4]
         return
     end
 
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/nlls_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/nlls_tests.jl
index dc3835b6..1a3bd892 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/nlls_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/nlls_tests.jl
@@ -12,8 +12,8 @@ SOLVERS_NAMES = ["$solver"
 # OOP MP-BVP
 f1(u, p, t) = [u[2], -u[1]]
 function bc1(sol, p, t)
-    solₜ₁ = sol[:, 1]
-    solₜ₂ = sol[:, end]
+    solₜ₁ = sol(0.0)
+    solₜ₂ = sol(100.0)
     return [solₜ₁[1], solₜ₂[1] - 1, solₜ₂[2] + 1.729109]
 end
 
@@ -25,8 +25,8 @@ function f1!(du, u, p, t)
 end
 
 function bc1!(resid, sol, p, t)
-    solₜ₁ = sol[:, 1]
-    solₜ₂ = sol[:, end]
+    solₜ₁ = sol(0.0)
+    solₜ₂ = sol(100.0)
     # We know that this overconstrained system has a solution
     resid[1] = solₜ₁[1]
     resid[2] = solₜ₂[1] - 1
@@ -111,8 +111,8 @@ function bc_b!(residual, y, p)
 end
 
 function bc!(residual, sol, p, t)
-    y1 = first(sol)
-    y2 = last(sol)
+    y1 = sol(0.0)
+    y2 = sol(0.5)
     R0_u = reshape(@view(y1[4:12]), 3, 3)
     RL_u = reshape(@view(y2[4:12]), 3, 3)
 
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/vectorofvector_initials_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/vectorofvector_initials_tests.jl
index b79d4461..1865f590 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/vectorofvector_initials_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/vectorofvector_initials_tests.jl
@@ -56,9 +56,9 @@
     # The BVP set up
     # This is not really kind of Two-Point BVP we support.
     function bc_po!(residual, u, p, t)
-        residual[1] = u[:, 1][1] - u[:, end][1]
-        residual[2] = u[:, 1][2] - u[:, end][2]
-        residual[3] = u[:, 1][3] - u[:, end][3]
+        residual[1] = u(0.0)[1] - u(T)[1]
+        residual[2] = u(0.0)[2] - u(T)[2]
+        residual[3] = u(0.0)[3] - u(T)[3]
     end
 
     #This is the part of the code that has problems
diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl
index 79331b7a..b0459673 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl
@@ -7,8 +7,8 @@
     end
 
     function bc!(residual, u, p, t)
-        residual[1] = u[:, 1][1] - 1.0
-        residual[2] = u[:, end][1]
+        residual[1] = u(0.0)[1] - 1.0
+        residual[2] = u(1.0)[1]
     end
 
     prob_func(prob, i, repeat) = remake(prob, p = [rand()])
diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
index ff7d4cf6..d4a9329a 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
@@ -39,10 +39,16 @@ end
 f2(u, p, t) = [u[2], -u[1]]
 
 function boundary!(residual, u, p, t)
+    residual[1] = u(0.0)[1] - 5
+    residual[2] = u(5.0)[1]
+end
+boundary(u, p, t) = [u(0.0)[1] - 5, u(5.0)[1]]
+
+function boundary_indexing!(residual, u, p, t)
     residual[1] = u[:, 1][1] - 5
     residual[2] = u[:, end][1]
 end
-boundary(u, p, t) = [u[:, 1][1] - 5, u[:, end][1]]
+boundary_indexing(u, p, t) = [u[:, 1][1] - 5, u[:, end][1]]
 
 function boundary_two_point_a!(resida, ua, p)
     resida[1] = ua[1] - 5
@@ -75,6 +81,8 @@ probArr = [BVProblem(odef1!, boundary!, u0, tspan, nlls = Val(false)),
     BVProblem(odef1, boundary, u0, tspan, nlls = Val(false)),
     BVProblem(odef2!, boundary!, u0, tspan, nlls = Val(false)),
     BVProblem(odef2, boundary, u0, tspan, nlls = Val(false)),
+    BVProblem(odef2!, boundary_indexing!, u0, tspan, nlls = Val(false)),
+    BVProblem(odef2, boundary_indexing, u0, tspan, nlls = Val(false)),
     TwoPointBVProblem(odef1!, (boundary_two_point_a!, boundary_two_point_b!),
         u0, tspan; bcresid_prototype, nlls = Val(false)),
     TwoPointBVProblem(odef1, (boundary_two_point_a, boundary_two_point_b),
@@ -96,7 +104,7 @@ end
 @testitem "Affineness" setup=[FIRKNestedConvergenceTests] begin
     using LinearAlgebra
 
-    @testset "Problem: $i" for i in (1, 2, 5, 6)
+    @testset "Problem: $i" for i in (1, 2, 7, 8)
         prob = probArr[i]
 
         @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5)
@@ -141,7 +149,7 @@ end
 @testitem "JET: Runtime Dispatches" setup=[FIRKNestedConvergenceTests] begin
     using JET
 
-    @testset "Problem: $i" for i in 1:8
+    @testset "Problem: $i" for i in 1:10
         prob = probArr[i]
         @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5)
             solver = lobattoIIIa_solver(Val(stage); nlsolve = NewtonRaphson(),
@@ -181,14 +189,14 @@ end
 @testitem "Convergence on Linear" setup=[FIRKNestedConvergenceTests] begin
     using LinearAlgebra, DiffEqDevTools
 
-    @testset "Problem: $i" for i in (3, 4, 7, 8)
+    @testset "Problem: $i" for i in (3, 4, 9, 10)
         prob = probArr[i]
 
         @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5)
             @time sim = test_convergence(
                 dts, prob, lobattoIIIa_solver(Val(stage); nested_nlsolve = nested);
                 abstol = 1e-8, reltol = 1e-8)
-            if (stage == 4 && ((i == 7) || (i == 8)))
+            if (stage == 4 && ((i == 9) || (i == 10)))
                 @test sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             elseif first(sim.errors[:final]) < 1e-12
                 @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
@@ -201,7 +209,7 @@ end
             @time sim = test_convergence(
                 dts, prob, lobattoIIIb_solver(Val(stage); nested_nlsolve = nested);
                 abstol = 1e-8, reltol = 1e-8)
-            if (stage == 4 && ((i == 7) || (i == 8)))
+            if (stage == 4 && ((i == 9) || (i == 10)))
                 @test sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
             elseif first(sim.errors[:final]) < 1e-12
                 @test_broken sim.𝒪est[:final]≈2 * stage - 2 atol=testTol
@@ -248,8 +256,8 @@ end
     end
 
     function bc_pendulum!(residual, u, p, t)
-        residual[1] = u[:, end ÷ 2][1] + π / 2 # the solution at the middle of the time span should be -pi/2
-        residual[2] = u[:, end][1] - π / 2 # the solution at the end of the time span should be pi/2
+        residual[1] = u(pi / 4)[1] + π / 2 # the solution at the middle of the time span should be -pi/2
+        residual[2] = u(pi / 2)[1] - π / 2 # the solution at the end of the time span should be pi/2
     end
 
     u0 = MVector{2}([pi / 2, pi / 2])
@@ -301,8 +309,8 @@ end
         du[2] = 1 / p * u[1]
     end
     function prob_bvp_linear_bc!(res, u, p, t)
-        res[1] = u[:, 1][1] - 1
-        res[2] = u[:, end][1]
+        res[1] = u(0.0)[1] - 1
+        res[2] = u(1.0)[1]
     end
 
     prob_bvp_linear_function = ODEFunction(
@@ -384,12 +392,12 @@ end
     end
 
     function swirling_flow_bc!(res, u, p, t)
-        res[1] = u[:, 1][1] + 1.0
-        res[2] = u[:, 1][3]
-        res[3] = u[:, 1][4]
-        res[4] = u[:, end][1] - 1.0
-        res[5] = u[:, end][3]
-        res[6] = u[:, end][4]
+        res[1] = u(0.0)[1] + 1.0
+        res[2] = u(0.0)[3]
+        res[3] = u(0.0)[4]
+        res[4] = u(1.0)[1] - 1.0
+        res[5] = u(1.0)[3]
+        res[6] = u(1.0)[4]
         return
     end
 
diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/nlls_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/nlls_tests.jl
index 7e7d1d10..e654ebba 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/nlls_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/nlls_tests.jl
@@ -15,8 +15,8 @@ nlsolve in ["NewtonRaphson", "GaussNewton", "TrustRegion"]]
 # OOP MP-BVP
 f1(u, p, t) = [u[2], -u[1]]
 function bc1(sol, p, t)
-    solₜ₁ = sol[:, 1]
-    solₜ₂ = sol[:, end]
+    solₜ₁ = sol(0.0)
+    solₜ₂ = sol(100.0)
     return [solₜ₁[1], solₜ₂[1] - 1, solₜ₂[2] + 1.729109]
 end
 
@@ -28,8 +28,8 @@ function f1!(du, u, p, t)
 end
 
 function bc1!(resid, sol, p, t)
-    solₜ₁ = sol[:, 1]
-    solₜ₂ = sol[:, end]
+    solₜ₁ = sol(0.0)
+    solₜ₂ = sol(100.0)
     # We know that this overconstrained system has a solution
     resid[1] = solₜ₁[1]
     resid[2] = solₜ₂[1] - 1
@@ -114,8 +114,8 @@ function bc_b!(residual, y, p)
 end
 
 function bc!(residual, sol, p, t)
-    y1 = first(sol)
-    y2 = last(sol)
+    y1 = sol(0.0)
+    y2 = sol(0.5)
     R0_u = reshape(@view(y1[4:12]), 3, 3)
     RL_u = reshape(@view(y2[4:12]), 3, 3)
 
@@ -184,7 +184,9 @@ end
     @testset "Problem: $i" for i in 1:2
         prob = UnderconstrainedProbArr[i]
         @testset "Solver: $name" for (name, solver) in zip(SOLVERS_NAMES, SOLVERS)
-            if name == "RadauIIa5 with GaussNewton"
+            if (i == 2) && ((name == "RadauIIa5 with GaussNewton") ||
+                (name == "RadauIIa5 with NewtonRaphson"))
+                # Actually have successful retcode
                 continue
             else
                 sol = solve(
diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/vectorofvector_initials_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/vectorofvector_initials_tests.jl
index 38e9f433..97f2bc1f 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/vectorofvector_initials_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/vectorofvector_initials_tests.jl
@@ -56,9 +56,9 @@
     # The BVP set up
     # This is not really kind of Two-Point BVP we support.
     function bc_po!(residual, u, p, t)
-        residual[1] = u[:, 1][1] - u[:, end][1]
-        residual[2] = u[:, 1][2] - u[:, end][2]
-        residual[3] = u[:, 1][3] - u[:, end][3]
+        residual[1] = u(0.0)[1] - u(T)[1]
+        residual[2] = u(0.0)[2] - u(T)[2]
+        residual[3] = u(0.0)[3] - u(T)[3]
     end
 
     bvp1 = BVProblem(TC!, bc_po!, sol.u, tspan)
diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
index 942e5b80..f271c207 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
@@ -15,7 +15,7 @@ import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorit
                                 recursive_flatten, recursive_flatten!, recursive_unflatten!,
                                 __concrete_nonlinearsolve_algorithm, diff!,
                                 __FastShortcutBVPCompatibleNonlinearPolyalg,
-                                __FastShortcutBVPCompatibleNLLSPolyalg,
+                                __FastShortcutBVPCompatibleNLLSPolyalg, __restructure_sol,
                                 concrete_jacobian_algorithm, eval_bc_residual,
                                 eval_bc_residual!, get_tmp, __maybe_matmul!,
                                 __append_similar!, __extract_problem_details,
diff --git a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
index 00f8db92..234e672b 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
@@ -117,8 +117,9 @@ end
 
 # Intermidiate solution for evaluating boundry conditions
 # basically simplified version of the interpolation for MIRK
-function (s::EvalSol)(tval::Number)
-    (; t, u, alg, k_discrete) = s
+function (s::EvalSol{C})(tval::Number) where {C <: MIRKCache}
+    (; t, u, cache) = s
+    (; alg, k_discrete) = cache
     stage = alg_stage(alg)
     # Quick handle for the case where tval is at the boundary
     (tval == t[1]) && return first(u)
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index 27ccf157..b69a8450 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -210,7 +210,7 @@ function __construct_nlproblem(
         cache::MIRKCache{iip}, y::AbstractVector, y₀::AbstractVectorOfArray) where {iip}
     pt = cache.problem_type
 
-    eval_sol = EvalSol{iip}(y₀.u, cache.mesh, cache.alg, cache.k_discrete)
+    eval_sol = EvalSol(__restructure_sol(y₀.u, cache.in_size), cache.mesh, cache)
 
     loss_bc = if iip
         @closure (du, u, p) -> __mirk_loss_bc!(
@@ -243,8 +243,8 @@ end
     y_ = recursive_unflatten!(y, u)
     resids = [get_tmp(r, u) for r in residual]
     Φ!(resids[2:end], cache, y_, u, p)
-    EvalSol.u[1:end] .= y_
-    EvalSol.k_discrete[1:end] .= cache.k_discrete
+    EvalSol.u[1:end] .= __restructure_sol(y_, cache.in_size)
+    EvalSol.cache.k_discrete[1:end] .= cache.k_discrete
     eval_bc_residual!(resids[1], pt, bc!, EvalSol, p, mesh)
     recursive_flatten!(resid, resids)
     return nothing
@@ -267,14 +267,14 @@ end
         u, p, y, pt::StandardBVProblem, bc::BC, mesh, cache, EvalSol) where {BC}
     y_ = recursive_unflatten!(y, u)
     resid_co = Φ(cache, y_, u, p)
-    EvalSol.u[1:end] .= y_
-    EvalSol.k_discrete[1:end] .= cache.k_discrete
+    EvalSol.u[1:end] .= __restructure_sol(y_, cache.in_size)
+    EvalSol.cache.k_discrete[1:end] .= cache.k_discrete
     resid_bc = eval_bc_residual(pt, bc, EvalSol, p, mesh)
     return vcat(resid_bc, mapreduce(vec, vcat, resid_co))
 end
 
 @views function __mirk_loss(u, p, y, pt::TwoPointBVProblem, bc::Tuple{BC1, BC2},
-        mesh, cache, EvalSol) where {BC1, BC2}
+        mesh, cache, _) where {BC1, BC2}
     y_ = recursive_unflatten!(y, u)
     resid_co = Φ(cache, y_, u, p)
     soly_ = VectorOfArray(y_)
@@ -285,14 +285,14 @@ end
 @views function __mirk_loss_bc!(
         resid, u, p, pt, bc!::BC, y, mesh, cache::MIRKCache) where {BC}
     y_ = recursive_unflatten!(y, u)
-    soly_ = EvalSol{true}(y_, mesh, cache.alg, cache.k_discrete)
+    soly_ = EvalSol(__restructure_sol(y_, cache.in_size), mesh, cache)
     eval_bc_residual!(resid, pt, bc!, soly_, p, mesh)
     return nothing
 end
 
 @views function __mirk_loss_bc(u, p, pt, bc!::BC, y, mesh, cache::MIRKCache) where {BC}
     y_ = recursive_unflatten!(y, u)
-    soly_ = EvalSol{false}(y_, mesh, cache.alg, cache.k_discrete)
+    soly_ = EvalSol(__restructure_sol(y_, cache.in_size), mesh, cache)
     return eval_bc_residual(pt, bc!, soly_, p, mesh)
 end
 
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index ffef123f..281e14fa 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -27,6 +27,13 @@ function boundary!(residual, u, p, t)
 end
 boundary(u, p, t) = [u(0.0)[1] - 5, u(5.0)[1]]
 
+# Array indexing for boudnary conditions
+function boundary_indexing!(residual, u, p, t)
+    residual[1] = u[:, 1][1] - 5
+    residual[2] = u[:, end][1]
+end
+boundary_indexing(u, p, t) = [u[:, 1][1] - 5, u[:, end][1]]
+
 function boundary_two_point_a!(resida, ua, p)
     resida[1] = ua[1] - 5
 end
@@ -58,6 +65,8 @@ probArr = [BVProblem(odef1!, boundary!, u0, tspan, nlls = Val(false)),
     BVProblem(odef1, boundary, u0, tspan, nlls = Val(false)),
     BVProblem(odef2!, boundary!, u0, tspan, nlls = Val(false)),
     BVProblem(odef2, boundary, u0, tspan, nlls = Val(false)),
+    BVProblem(odef2!, boundary_indexing!, u0, tspan, nlls = Val(false)),
+    BVProblem(odef2, boundary_indexing, u0, tspan, nlls = Val(false)),
     TwoPointBVProblem(odef1!, (boundary_two_point_a!, boundary_two_point_b!),
         u0, tspan; bcresid_prototype, nlls = Val(false)),
     TwoPointBVProblem(odef1, (boundary_two_point_a, boundary_two_point_b),
@@ -78,7 +87,7 @@ end
 @testitem "Affineness" setup=[MIRKConvergenceTests] begin
     using LinearAlgebra
 
-    @testset "Problem: $i" for i in (1, 2, 5, 6)
+    @testset "Problem: $i" for i in (1, 2, 7, 8)
         prob = probArr[i]
         @testset "MIRK$order" for order in (2, 3, 4, 5, 6)
             sol = solve(prob, mirk_solver(Val(order)); dt = 0.2)
@@ -90,7 +99,7 @@ end
 @testitem "JET: Runtime Dispatches" setup=[MIRKConvergenceTests] begin
     using JET
 
-    @testset "Problem: $i" for i in 1:8
+    @testset "Problem: $i" for i in 1:10
         prob = probArr[i]
         @testset "MIRK$order" for order in (2, 3, 4, 5, 6)
             solver = mirk_solver(Val(order); nlsolve = NewtonRaphson(),
@@ -106,7 +115,7 @@ end
 @testitem "Convergence on Linear" setup=[MIRKConvergenceTests] begin
     using LinearAlgebra, DiffEqDevTools
 
-    @testset "Problem: $i" for i in (3, 4, 7, 8)
+    @testset "Problem: $i" for i in (3, 4, 5, 6, 9, 10)
         prob = probArr[i]
         @testset "MIRK$order" for (_, order) in enumerate((2, 3, 4, 5, 6))
             sim = test_convergence(
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
index 9dde9815..aec4d55f 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
@@ -16,17 +16,18 @@ import BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorit
                                 __concrete_nonlinearsolve_algorithm, diff!,
                                 __FastShortcutBVPCompatibleNonlinearPolyalg,
                                 __FastShortcutBVPCompatibleNLLSPolyalg, eval_bc_residual,
-                                eval_bc_residual!, get_tmp, __maybe_matmul!,
+                                eval_bc_residual!, get_tmp, __maybe_matmul!, EvalSol,
                                 __append_similar!, __extract_problem_details,
                                 __initial_guess, __maybe_allocate_diffcache,
-                                __get_bcresid_prototype, __similar, __vec, __vec_f,
-                                __vec_f!, __vec_bc, __vec_bc!, recursive_flatten_twopoint!,
-                                __internal_nlsolve_problem, __extract_mesh, __extract_u0,
-                                __has_initial_guess, __initial_guess_length,
-                                __initial_guess_on_mesh, __flatten_initial_guess,
-                                __build_solution, __Fix3, __sparse_jacobian_cache,
-                                __sparsity_detection_alg, _sparse_like, ColoredMatrix,
-                                __default_sparse_ad, __default_nonsparse_ad
+                                __restructure_sol, __get_bcresid_prototype, __similar,
+                                __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
+                                recursive_flatten_twopoint!, __internal_nlsolve_problem,
+                                __extract_mesh, __extract_u0, __has_initial_guess,
+                                __initial_guess_length, __initial_guess_on_mesh,
+                                __flatten_initial_guess, __build_solution, __Fix3,
+                                __sparse_jacobian_cache, __sparsity_detection_alg,
+                                _sparse_like, ColoredMatrix, __default_sparse_ad,
+                                __default_nonsparse_ad
 
 import ADTypes: AbstractADType
 import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
@@ -48,6 +49,7 @@ include("mirkn.jl")
 include("alg_utils.jl")
 include("collocation.jl")
 include("mirkn_tableaus.jl")
+include("interpolation.jl")
 
 function __solve(
         prob::AbstractBVProblem, alg::BoundaryValueDiffEqAlgorithm, args...; kwargs...)
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/interpolation.jl b/lib/BoundaryValueDiffEqMIRKN/src/interpolation.jl
new file mode 100644
index 00000000..ad63c44d
--- /dev/null
+++ b/lib/BoundaryValueDiffEqMIRKN/src/interpolation.jl
@@ -0,0 +1,15 @@
+# Intermidiate solution for evaluating boundry conditions
+# basically simplified version of the linear interpolation for MIRKN
+function (s::EvalSol{C})(tval::Number) where {C <: MIRKNCache}
+    (; t, u, cache) = s
+
+    # Quick handle for the case where tval is at the boundary
+    (tval == t[1]) && return first(u)
+    (tval == t[end]) && return last(u)
+    # Linear interpolation
+    i = interval(tval, t)
+    dt = t[i + 1] - t[i]
+    τ = (tval - t[i]) / dt
+    z = τ * u[i + 1] + (1 - τ) * u[i]
+    return z
+end
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
index e9a753ef..3909b0f2 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
@@ -25,8 +25,8 @@ end
 
 Base.eltype(::MIRKNCache{iip, T}) where {iip, T} = T
 
-function SciMLBase.__init(
-        prob::SecondOrderBVProblem, alg::AbstractMIRKN; dt = 0.0, kwargs...)
+function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN;
+        dt = 0.0, adaptive = false, kwargs...)
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
     iip = isinplace(prob)
     t₀, t₁ = prob.tspan
@@ -160,8 +160,12 @@ end
         bc::BC, residual, mesh, cache::MIRKNCache) where {BC}
     y_ = recursive_unflatten!(y, u)
     resids = [get_tmp(r, u) for r in residual]
-    eval_bc_residual!(resids[1:2], pt, bc, y_, p, mesh)
     Φ!(resids[3:end], cache, y_, u, p)
+    soly_ = EvalSol(
+        __restructure_sol(y_[1:length(cache.mesh)], cache.in_size), cache.mesh, cache)
+    dsoly_ = EvalSol(__restructure_sol(y_[(length(cache.mesh) + 1):end], cache.in_size),
+        cache.mesh, cache)
+    eval_bc_residual!(resids[1:2], pt, bc, soly_, dsoly_, p, mesh)
     recursive_flatten!(resid, resids)
     return nothing
 end
@@ -169,8 +173,12 @@ end
 @views function __mirkn_loss(u, p, y, pt::StandardSecondOrderBVProblem,
         bc::BC, mesh, cache::MIRKNCache) where {BC}
     y_ = recursive_unflatten!(y, u)
-    resid_bc = eval_bc_residual(pt, bc, y_, p, mesh)
     resid_co = Φ(cache, y_, u, p)
+    soly_ = EvalSol(
+        __restructure_sol(y_[1:length(cache.mesh)], cache.in_size), cache.mesh, cache)
+    dsoly_ = EvalSol(__restructure_sol(y_[(length(cache.mesh) + 1):end], cache.in_size),
+        cache.mesh, cache)
+    resid_bc = eval_bc_residual(pt, bc, soly_, dsoly_, p, mesh)
     return vcat(resid_bc, mapreduce(vec, vcat, resid_co))
 end
 
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/utils.jl b/lib/BoundaryValueDiffEqMIRKN/src/utils.jl
index 630df39d..809d29b1 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/utils.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/utils.jl
@@ -5,9 +5,8 @@ function general_eval_bc_residual(
     res_bc = bc(y[(L + 1):end], y[1:L], p, mesh)
     return res_bc
 end
-function eval_bc_residual(::StandardSecondOrderBVProblem, bc::BC, y, p, mesh) where {BC}
-    L = length(mesh)
-    res_bc = bc(y[(L + 1):end], y[1:L], p, mesh)
+function eval_bc_residual(::StandardSecondOrderBVProblem, bc::BC, y, dy, p, mesh) where {BC}
+    res_bc = bc(dy, y, p, mesh)
     return res_bc
 end
 function eval_bc_residual(::TwoPointSecondOrderBVProblem, bc::BC, sol, p, mesh) where {BC}
@@ -20,25 +19,16 @@ function eval_bc_residual(::TwoPointSecondOrderBVProblem, bc::BC, sol, p, mesh)
     return vcat(bc[1](dua, ua, p), bc[2](dub, ub, p))
 end
 
-function eval_bc_residual!(resid, ::StandardBVProblem, bc!::BC, sol, p, t) where {BC}
-    bc!(resid, sol, p, t)
-end
-function eval_bc_residual!(
-        resid, ::StandardSecondOrderBVProblem, bc!::BC, dsol, sol, p, t) where {BC}
-    bc!(resid, dsol, sol, p, t)
-end
-
 function general_eval_bc_residual!(
         resid, ::StandardSecondOrderBVProblem, bc!::BC, y, p, mesh) where {BC}
     L = length(mesh)
     bc!(resid, y[(L + 1):end], y[1:L], p, mesh)
 end
 function eval_bc_residual!(
-        resid, ::StandardSecondOrderBVProblem, bc!::BC, y, p, mesh) where {BC}
-    M = length(y[1])
-    L = length(mesh)
+        resid, ::StandardSecondOrderBVProblem, bc!::BC, y, dy, p, mesh) where {BC}
+    M = length(resid[1])
     res_bc = vcat(resid[1], resid[2])
-    bc!(res_bc, y[(L + 1):end], y[1:L], p, mesh)
+    bc!(res_bc, dy, y, p, mesh)
     copyto!(resid[1], res_bc[1:M])
     copyto!(resid[2], res_bc[(M + 1):end])
 end
diff --git a/lib/BoundaryValueDiffEqMIRKN/test/mirkn_basic_tests.jl b/lib/BoundaryValueDiffEqMIRKN/test/mirkn_basic_tests.jl
index 8f309323..fdb83c1b 100644
--- a/lib/BoundaryValueDiffEqMIRKN/test/mirkn_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/test/mirkn_basic_tests.jl
@@ -18,17 +18,32 @@
     end
 
     function bc!(res, du, u, p, t)
-        res[1] = u[1][1]
-        res[2] = u[end][1]
-        res[3] = u[1][3] + 1
-        res[4] = u[end][3] - 1
-        res[5] = du[1][1]
-        res[6] = du[end][1]
+        res[1] = u(0.0)[1]
+        res[2] = u(1.0)[1]
+        res[3] = u(0.0)[3] + 1
+        res[4] = u(1.0)[3] - 1
+        res[5] = du(0.0)[1]
+        res[6] = du(1.0)[1]
     end
 
     function bc(du, u, p, t)
-        return [u[1][1], u[end][1], u[1][3] + 1, u[end][3] - 1, du[1][1], du[end][1]]
+        return [u(0.0)[1], u(1.0)[1], u(0.0)[3] + 1, u(1.0)[3] - 1, du(0.0)[1], du(1.0)[1]]
     end
+
+    function bc_indexing!(res, du, u, p, t)
+        res[1] = u[:, 1][1]
+        res[2] = u[:, end][1]
+        res[3] = u[:, 1][3] + 1
+        res[4] = u[:, end][3] - 1
+        res[5] = du[:, 1][1]
+        res[6] = du[:, end][1]
+    end
+
+    function bc_indexing(du, u, p, t)
+        return [u[:, 1][1], u[:, end][1], u[:, 1][3] + 1,
+            u[:, end][3] - 1, du[:, 1][1], du[:, end][1]]
+    end
+
     function bca!(resa, du, u, p)
         resa[1] = u[1]
         resa[2] = u[3] + 1
@@ -52,13 +67,15 @@
 
     probArr = [SecondOrderBVProblem(test!, bc!, u0, tspan),
         SecondOrderBVProblem(test, bc, u0, tspan),
+        SecondOrderBVProblem(test!, bc_indexing!, u0, tspan),
+        SecondOrderBVProblem(test, bc_indexing, u0, tspan),
         TwoPointSecondOrderBVProblem(
             test!, (bca!, bcb!), u0, tspan, bcresid_prototype = (zeros(3), zeros(3))),
         TwoPointSecondOrderBVProblem(
             test, (bca, bcb), u0, tspan, bcresid_prototype = (zeros(3), zeros(3)))]
 
     @testset "MIRKN$order" for order in (4, 6)
-        @testset "Problem $i" for i in 1:4
+        @testset "Problem $i" for i in 1:6
             sol = solve(probArr[i], mirkn_solver(Val(order)); dt = 0.01)
             @test SciMLBase.successful_retcode(sol)
         end
@@ -80,11 +97,18 @@ end
         return u[1]
     end
     function bc!(res, du, u, p, t)
-        res[1] = u[1][1] - 1
-        res[2] = u[end][1]
+        res[1] = u(0.0)[1] - 1
+        res[2] = u(1.0)[1]
     end
     function bc(du, u, p, t)
-        return [u[1][1] - 1, u[end][1]]
+        return [u(0.0)[1] - 1, u(1.0)[1]]
+    end
+    function bc_indexing!(res, du, u, p, t)
+        res[1] = u[:, 1][1] - 1
+        res[2] = u[:, end][1]
+    end
+    function bc_indexing(du, u, p, t)
+        return [u[:, 1][1] - 1, u[:, end][1]]
     end
     function bc_a!(res, du, u, p)
         res[1] = u[1] - 1
@@ -98,35 +122,26 @@ end
     function bc_b(du, u, p)
         return [u[1]]
     end
+    analytical_solution = (u0, p, t) -> [
+        (exp(-t) - exp(t - 2)) / (1 - exp(-2)), (-exp(-t) - exp(t - 2)) / (1 - exp(-2))]
     u0 = [1.0]
     tspan = (0.0, 1.0)
     testTol = 0.2
-    bvpf1 = DynamicalBVPFunction(f!,
-        bc!,
-        analytic = (u0, p, t) -> [
-            (exp(-t) - exp(t - 2)) / (1 - exp(-2)), (-exp(-t) - exp(t - 2)) / (1 - exp(-2))])
-    bvpf2 = DynamicalBVPFunction(f,
-        bc,
-        analytic = (u0, p, t) -> [
-            (exp(-t) - exp(t - 2)) / (1 - exp(-2)), (-exp(-t) - exp(t - 2)) / (1 - exp(-2))])
-    bvpf3 = DynamicalBVPFunction(f!,
-        (bc_a!, bc_b!),
-        analytic = (u0, p, t) -> [
-            (exp(-t) - exp(t - 2)) / (1 - exp(-2)), (-exp(-t) - exp(t - 2)) / (1 - exp(-2))],
-        bcresid_prototype = (zeros(1), zeros(1)),
-        twopoint = Val(true))
-    bvpf4 = DynamicalBVPFunction(f,
-        (bc_a, bc_b),
-        analytic = (u0, p, t) -> [
-            (exp(-t) - exp(t - 2)) / (1 - exp(-2)), (-exp(-t) - exp(t - 2)) / (1 - exp(-2))],
-        bcresid_prototype = (zeros(1), zeros(1)),
-        twopoint = Val(true))
+    bvpf1 = DynamicalBVPFunction(f!, bc!, analytic = analytical_solution)
+    bvpf2 = DynamicalBVPFunction(f, bc, analytic = analytical_solution)
+    bvpf3 = DynamicalBVPFunction(f!, bc_indexing!, analytic = analytical_solution)
+    bvpf4 = DynamicalBVPFunction(f, bc_indexing, analytic = analytical_solution)
+    bvpf5 = DynamicalBVPFunction(f!, (bc_a!, bc_b!), analytic = analytical_solution,
+        bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
+    bvpf6 = DynamicalBVPFunction(f, (bc_a, bc_b), analytic = analytical_solution,
+        bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
     probArr = [
         SecondOrderBVProblem(bvpf1, u0, tspan), SecondOrderBVProblem(bvpf2, u0, tspan),
-        TwoPointSecondOrderBVProblem(bvpf3, u0, tspan),
-        TwoPointSecondOrderBVProblem(bvpf4, u0, tspan)]
+        SecondOrderBVProblem(bvpf3, u0, tspan), SecondOrderBVProblem(bvpf4, u0, tspan),
+        TwoPointSecondOrderBVProblem(bvpf5, u0, tspan),
+        TwoPointSecondOrderBVProblem(bvpf6, u0, tspan)]
     dts = 1 .// 2 .^ (3:-1:1)
-    @testset "Problem: $i" for i in (1, 2, 3, 4)
+    @testset "Problem: $i" for i in (1, 2, 3, 4, 5, 6)
         prob = probArr[i]
         @testset "MIRKN$order" for order in (4, 6)
             sim = test_convergence(
diff --git a/test/misc/bigfloat_test.jl b/test/misc/bigfloat_test.jl
index c202bf78..a368dd83 100644
--- a/test/misc/bigfloat_test.jl
+++ b/test/misc/bigfloat_test.jl
@@ -8,8 +8,8 @@
         du[2] = -9.81 * sin(θ)
     end
     function bc!(residual, u, p, t)
-        residual[1] = u[:, end ÷ 2][1] + big(pi / 2)
-        residual[2] = u[:, end][1] - big(pi / 2)
+        residual[1] = u(pi / 4)[1] + big(pi / 2)
+        residual[2] = u(pi / 2)[1] - big(pi / 2)
     end
     u0 = BigFloat.([pi / 2, pi / 2])
     prob = BVProblem(simplependulum!, bc!, u0, tspan)
diff --git a/test/misc/manifolds_tests.jl b/test/misc/manifolds_tests.jl
index 42a5b253..d47a5340 100644
--- a/test/misc/manifolds_tests.jl
+++ b/test/misc/manifolds_tests.jl
@@ -27,9 +27,9 @@
 
     function bc1!(residual, u, params, t)
         M, i, a1, a2 = params
-        mid = div(length(u[:, 1]), 2)
-        residual[1:mid] = u[:, 1][1:mid] - a1
-        residual[(mid + 1):end] = u[:, end][1:mid] - a2
+        mid = div(length(u(0.0)), 2)
+        residual[1:mid] = u(0.0)[1:mid] - a1
+        residual[(mid + 1):end] = u(1.0)[1:mid] - a2
         return
     end
 
@@ -67,7 +67,7 @@
             if alg isa Shooting || alg isa MultipleShooting
                 sol = solve(bvp, alg)
             else
-                sol = solve(bvp, alg; dt)
+                sol = solve(bvp, alg; dt, abstol = 1e-8)
             end
             @test SciMLBase.successful_retcode(sol)
             resid = zeros(4)
diff --git a/test/misc/non_vector_input_tests.jl b/test/misc/non_vector_input_tests.jl
index 5c6cbf1f..95413aa9 100644
--- a/test/misc/non_vector_input_tests.jl
+++ b/test/misc/non_vector_input_tests.jl
@@ -16,8 +16,8 @@
     end
 
     function boundary!(residual, u, p, t)
-        residual[1, 1] = u.u[1][1, 1] - 5
-        residual[1, 2] = u.u[end][1, 1]
+        residual[1, 1] = u(0.0)[1, 1] - 5
+        residual[1, 2] = u(5.0)[1, 1]
     end
 
     function boundary_a!(resida, ua, p)
@@ -28,7 +28,7 @@
     end
 
     function boundary(u, p, t)
-        return [u.u[1][1, 1] - 5 u.u[end][1, 1]]
+        return [u(0.0)[1, 1] - 5 u(5.0)[1, 1]]
     end
 
     boundary_a = (ua, p) -> [ua[1, 1] - 5]
diff --git a/test/misc/type_stability_tests.jl b/test/misc/type_stability_tests.jl
index 262b296f..19a82793 100644
--- a/test/misc/type_stability_tests.jl
+++ b/test/misc/type_stability_tests.jl
@@ -7,10 +7,10 @@
         du[2] = p[3] * u[1] * u[2] - p[4] * u[2]
     end
 
-    bc(sol, p, t) = [sol[:, 1][1] - 1, sol[:, end][2] - 2]
+    bc(sol, p, t) = [sol(0.0)[1] - 1, sol(1.0)[2] - 2]
     function bc!(res, sol, p, t)
-        res[1] = sol[:, 1][1] - 1
-        res[2] = sol[:, end][2] - 2
+        res[1] = sol(0.0)[1] - 1
+        res[2] = sol(1.0)[2] - 2
     end
     twobc_a(ua, p) = [ua[1] - 1]
     twobc_b(ub, p) = [ub[2] - 2]