From 9b38ab170221ffdbfdbe3129119167049e6ebf39 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?=
 <sebastian.mc95@proton.me>
Date: Tue, 24 Mar 2026 02:50:26 +0200
Subject: [PATCH 1/2] test(bvp): remove direct `EvalSol` comparison

The discrepancy between the interpolation
methods is fixed by
https://github.com/SciML/BoundaryValueDiffEq.jl/pull/452
---
 lib/ModelingToolkitBase/test/bvproblem.jl | 17 +++++------------
 1 file changed, 5 insertions(+), 12 deletions(-)

diff --git a/lib/ModelingToolkitBase/test/bvproblem.jl b/lib/ModelingToolkitBase/test/bvproblem.jl
index 01756e0a2b..e320f6da3c 100644
--- a/lib/ModelingToolkitBase/test/bvproblem.jl
+++ b/lib/ModelingToolkitBase/test/bvproblem.jl
@@ -149,24 +149,17 @@ end
 end
 
 function test_solvers(
-        solvers, prob, u0map, constraints, equations = []; dt = 0.005, atol = 1.0e-2
+        solvers, prob, u0map, constraints, equations = []; dt = 0.005, atol = 5.0e-6
     )
     for solver in solvers
         println("Solver: $solver")
-        sol = solve(prob, solver(), dt = dt, abstol = 1.0e-8, reltol = 1.0e-8)
+        sol = solve(prob, solver(), dt = dt)
         @test SciMLBase.successful_retcode(sol.retcode)
         p = prob.p
         t = sol.t
         bc = prob.f.bc
         ns = length(prob.u0)
 
-        cache = init(prob, solver(), dt = dt, abstol = 1.0e-8, reltol = 1.0e-8)
-        eval_sol = BoundaryValueDiffEqMIRK.EvalSol(
-            BoundaryValueDiffEqMIRK.__restructure_sol(
-                sol.u, cache.in_size
-            ), cache.mesh, cache
-        )
-        @test bc(zeros(ns), eval_sol, p, t) ≈ zeros(ns)
         if isinplace(prob.f)
             resid = zeros(ns)
             bc(resid, sol, p, t)
@@ -211,7 +204,7 @@ end
     bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(
         lksys, u0map, tspan; guesses = guess, cse = false
     )
-    test_solvers(solvers, bvp, u0map, constr; dt = 1.0e-4)
+    test_solvers(solvers, bvp, u0map, constr; dt = 1.0e-2)
 
     # Testing that more complicated constraints give correct solutions.
     constr = [EvalAt(0.2)(y) + EvalAt(0.8)(x) ~ 3.0, EvalAt(0.3)(y) ~ 2.0]
@@ -219,14 +212,14 @@ end
     bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(
         lksys, u0map, tspan; guesses = guess
     )
-    test_solvers(solvers, bvp, u0map, constr; dt = 1.0e-4)
+    test_solvers(solvers, bvp, u0map, constr; dt = 1.0e-2)
 
     constr = [α * β - EvalAt(0.6)(x) ~ 0.0, EvalAt(0.2)(y) ~ 3.0]
     @mtkcompile lksys = System(eqs, t; constraints = constr)
     bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(
         lksys, u0map, tspan; guesses = guess
     )
-    test_solvers(solvers, bvp, u0map, constr; dt = 1.0e-4)
+    test_solvers(solvers, bvp, u0map, constr; dt = 1.0e-2)
 end
 
 # Cartesian pendulum from the docs.

From d3b929f10af5ce078974ea4ffa8346d6c7804492 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?=
 <sebastian.mc95@proton.me>
Date: Tue, 24 Mar 2026 02:52:27 +0200
Subject: [PATCH 2/2] test(bvp): compare the solution at stricter tolerances

for the BVP case the solution is more precise at default tolerances
---
 lib/ModelingToolkitBase/test/bvproblem.jl | 16 ++++++----------
 1 file changed, 6 insertions(+), 10 deletions(-)

diff --git a/lib/ModelingToolkitBase/test/bvproblem.jl b/lib/ModelingToolkitBase/test/bvproblem.jl
index e320f6da3c..df37b7fff2 100644
--- a/lib/ModelingToolkitBase/test/bvproblem.jl
+++ b/lib/ModelingToolkitBase/test/bvproblem.jl
@@ -67,25 +67,21 @@ end
     tspan = (0.0, 6.0)
 
     op = ODEProblem(pend, [u0map; parammap], tspan)
-    osol = solve(op, Vern9())
+    osol = solve(op, Vern9(), abstol=1e-10, reltol=1e-10)
 
-    bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(
-        pend, [u0map; parammap], tspan
-    )
+    bvp = BVProblem(pend, [u0map; parammap], tspan)
     for solver in solvers
-        sol = solve(bvp, solver(), dt = 0.01)
-        @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
+        sol = solve(bvp, solver(), dt = 1e-2)
+        @test isapprox(sol.u[end], osol.u[end])
         @test sol.u[1] == [π / 2, π / 2]
     end
 
     # Test out-of-place
-    bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(
-        pend, [u0map; parammap], tspan
-    )
+    bvp2 = BVProblem{false, SciMLBase.FullSpecialize}(pend, [u0map; parammap], tspan)
 
     for solver in solvers
         sol = solve(bvp2, solver(), dt = 0.01)
-        @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
+        @test isapprox(sol.u[end], osol.u[end])
         @test sol.u[1] == [π / 2, π / 2]
     end
 end