diff --git a/src/components.jl b/src/components.jl
index 3d6c2676..8f9b67d3 100644
--- a/src/components.jl
+++ b/src/components.jl
@@ -540,6 +540,7 @@ See also [`BodyCylinder`](@ref) and [`BodyBox`](@ref) for body components with p
         frame_cm = Frame()
     end
 
+    # NOTE: these parameters should be defined before the `systems` block above, but due to bugs in MTK/JSC with higher-order array parameters we cannot do that. We still define the parameters so that they are available to make animations
     @variables r_0(t)[1:3]=r_0 [
         state_priority = 2,
         description = "Position vector from origin of world frame to origin of frame_a",
diff --git a/test/runtests.jl b/test/runtests.jl
index 846f072e..0124a47c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,8 +4,8 @@ using Test
 using JuliaSimCompiler
 using OrdinaryDiffEq
 using LinearAlgebra
-t = Multibody.t
-D = Differential(t)
+isdefined(Main, :t) || (t = Multibody.t)
+isdefined(Main, :D) || (D = Differential(t))
 doplot() = false
 world = Multibody.world
 
@@ -1124,6 +1124,7 @@ using LinearAlgebra
 
     prob = ODEProblem(ssys, [model.joint.phi => 0], (0, 1))
     sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
+    @test SciMLBase.successful_retcode(sol)
     # first(render(model, sol, 0, x=2.5, y=1.5, z=2.5, show_axis=true))
     # @test sol(10, idxs=model.body.body.m) ≈ 226.27 rtol=1e-3 # Values from open modelica
     # @test sol(10, idxs=model.body.body.I_11) ≈ 245.28 rtol=1e-3
@@ -1148,8 +1149,6 @@ end
 
 using LinearAlgebra, ModelingToolkit, Multibody, JuliaSimCompiler, OrdinaryDiffEq
 using Multibody.Rotations: RotXYZ
-t = Multibody.t
-D = Multibody.D
 world = Multibody.world
 
 @named joint = Multibody.Spherical(isroot=false, state=false, quat=false)
diff --git a/test/test_quaternions.jl b/test/test_quaternions.jl
index 28c87b64..103d5ccb 100644
--- a/test/test_quaternions.jl
+++ b/test/test_quaternions.jl
@@ -198,7 +198,6 @@ end
 
 # @testset "Spherical joint with quaternion state" begin
     using LinearAlgebra, ModelingToolkit, Multibody, JuliaSimCompiler
-    t = Multibody.t
     world = Multibody.world
 
 
@@ -220,7 +219,7 @@ end
     ssys = structural_simplify(irsys)
 
 
-    D = Differential(t)
+    isdefined(Main, :D) || (D = Differential(t))
     # q0 = randn(4); q0 ./= norm(q0)
     # q0 = [1,0,0,0]
     prob = ODEProblem(ssys, [