Differences between the interpolation between the solution and the intermediarry `EvalSol`

**Describe the bug 🐞** Boundary conditions that were satisfied to machine precision during the solve can show non-trivial residuals when re-evaluated on the returned solution object via `sol(t)`. Relevant source locations (BoundaryValueDiffEqMIRK): - `EvalSol` interpolation: `src/interpolation.jl:205-231` (uses `evalsol_interp_weights`, `k_discrete` only) - `MIRKInterpolation`: `src/interpolation.jl:93-150` (uses `interp_weights`, `k_discrete` + `k_interp`) - `EvalSol` struct: `BoundaryValueDiffEqCore/src/solution_utils.jl:2-6` **Expected behavior** A clear and concise description of what you expected to happen. **Minimal Reproducible Example 👇** Without MRE, we would only be able to help you to a limited extent, and attention to the issue would be limited. to know more about MRE refer to [wikipedia](https://en.wikipedia.org/wiki/Minimal_reproducible_example) and [stackoverflow](https://stackoverflow.com/help/minimal-reproducible-example). ```julia function lotkavolterra!(du, u, p, t) α, β, δ, γ = p du[1] = α * u[1] - β * u[1] * u[2] du[2] = -γ * u[2] + δ * u[1] * u[2] end # Boundary conditions: x(0.6) = 3.5, x(0.3) = 7.0 # These are interior-point conditions — they require interpolation. function bc2!(resid, u, p, t) resid[1] = u(0.6)[1] - 3.5 resid[2] = u(0.3)[1] - 7.0 end u0_manual = [4.0, 2.0] tspan_manual = (0.0, 1.0) p_manual = [1.5, 1.0, 1.0, 3.0] prob2 = BVProblem(lotkavolterra!, bc2!, u0_manual, tspan_manual, p_manual) dt2 = 0.005 sol2 = solve(prob2, MIRK4(), dt = dt2, abstol = 1e-8, reltol = 1e-8) println("Solver retcode: ", sol2.retcode) # --- Reconstruct EvalSol from solver internals --- cache2 = init(prob2, MIRK4(), dt = dt2, abstol = 1e-8, reltol = 1e-8) eval_sol2 = BoundaryValueDiffEqMIRK.EvalSol( BoundaryValueDiffEqMIRK.__restructure_sol(sol2.u, cache2.in_size), cache2.mesh, cache2 ) # --- Compare BC residuals --- resid_evalsol2 = zeros(2) bc2!(resid_evalsol2, eval_sol2, p_manual, sol2.t) resid_sol2 = zeros(2) bc2!(resid_sol2, sol2, p_manual, sol2.t) println("\n=== BVP Interpolation Discrepancy (no MTK) ===") println("Constraints: u[1](0.6) = 3.5, u[1](0.3) = 7.0") println() println("BC residuals using EvalSol (what the solver optimized against):") println(" |resid| = ", norm(resid_evalsol2)) println(" resid = ", resid_evalsol2) println() println("BC residuals using sol (what the user gets back):") println(" |resid| = ", norm(resid_sol2)) println(" resid = ", resid_sol2) println() println("=== Interpolation values at constraint points ===") for t_eval in [0.3, 0.6] println(" t = $t_eval:") println(" EvalSol($t_eval) = ", eval_sol2(t_eval)) println(" sol($t_eval) = ", sol2(t_eval)) println(" difference = ", eval_sol2(t_eval) .- sol2(t_eval)) end println() println("=== Interpolation difference across the domain ===") println(" t | ‖EvalSol(t) - sol(t)‖") println(" --------|----------------------") for t_eval in 0.1:0.1:0.9 diff = norm(eval_sol2(t_eval) .- sol2(t_eval)) println(" $(lpad(t_eval, 7)) | $diff") end ``` **Error & Stacktrace ⚠️** ``` Constraints: u[1](0.6) = 3.5, u[1](0.3) = 7.0 BC residuals using EvalSol (what the solver optimized against): |resid| = 0.0 resid = [0.0, 0.0] BC residuals using sol (what the user gets back): |resid| = 0.06646094374232514 resid = [-0.05862682879409942, -0.0313041848427309] === Interpolation values at constraint points === t = 0.3: EvalSol(0.3) = [7.0, 2.3725883756653507] sol(0.3) = [6.968695815157269, 2.4203299566370995] difference = [0.0313041848427309, -0.04774158097174874] t = 0.6: EvalSol(0.6) = [3.5, 4.8726232553145685] sol(0.6) = [3.4413731712059006, 4.884102541830545] difference = [0.05862682879409942, -0.011479286515976739] === Interpolation difference across the domain === t | ‖EvalSol(t) - sol(t)‖ --------|---------------------- 0.1 | 0.027232864467045273 0.2 | 0.03408968174242 0.3 | 0.05708949590204756 0.4 | 0.07862356204274779 0.5 | 0.07849837772348761 0.6 | 0.0597400960274465 0.7 | 0.043310821446036214 0.8 | 0.03751369110888054 0.9 | 0.035104877243346905 ``` **Environment (please complete the following information):** - Output of `using Pkg; Pkg.status()` ```julia ``` - Output of `using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)` ```julia ``` - Output of `versioninfo()` ```julia ``` **Additional context** Add any other context about the problem here.