diff --git a/benchmarks/ModelingToolkit/Multibody_Robot.jmd b/benchmarks/ModelingToolkit/Multibody_Robot.jmd
index b4de60dd9..13954afc8 100644
--- a/benchmarks/ModelingToolkit/Multibody_Robot.jmd
+++ b/benchmarks/ModelingToolkit/Multibody_Robot.jmd
@@ -74,6 +74,7 @@ f
 using OMJulia
 mod = OMJulia.OMCSession();
 OMJulia.sendExpression(mod, "getVersion()")
+OMJulia.sendExpression(mod, "installPackage(Modelica)")
 
 @show "Start OpenModelica Timings"
 
diff --git a/benchmarks/ModelingToolkit/Project.toml b/benchmarks/ModelingToolkit/Project.toml
index 0b71b5ef5..27361e1c8 100644
--- a/benchmarks/ModelingToolkit/Project.toml
+++ b/benchmarks/ModelingToolkit/Project.toml
@@ -6,6 +6,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d"
 JuliaSimCompilerRuntime = "9cbdfd5a-25c0-4dde-8b55-206f91b28bd9"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+Multibody = "e1cad5d1-98ef-44f9-a79a-9ca4547f95b9"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 OMJulia = "0f4fe800-344e-11e9-2949-fb537ad918e1"
@@ -24,6 +25,7 @@ CSV = "0.10"
 JuliaSimCompiler = "0.1.9"
 JuliaSimCompilerRuntime = "1.0.0"
 LinearSolve = "2.30.0"
+Multibody = "0.1.0"
 ModelingToolkit = "9.19.0"
 ModelingToolkitStandardLibrary = "2"
 OMJulia = "0.3.1"
@@ -31,4 +33,4 @@ OrdinaryDiffEq = "6.84.0"
 Polynomials = "4.0.8"
 PreferenceTools = "0.1.2"
 SciMLBenchmarks = "0.1.3"
-Symbolics = "5.30.1"
+Symbolics = "5.30.4"
diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 2808becd5..d4d78a594 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -93,12 +93,12 @@ end
 N = [5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 20000];
 
 # max size we test per method
-max_sizes = [4_000, 8_000, 20_000, 20_000, 20_000, 20_000, 20_000, 9000, 20_000];
+max_sizes = [4_000, 8_000, 20_000, 20_000, 20_000, 20_000, 20_000, 9000];
 
 # NaN-initialize so Makie will ignore incomplete
 ss_times = fill(NaN, length(N), 3);
 times = fill((NaN,NaN,NaN), length(N), length(max_sizes) - 1);
-total_times = fill(NaN, length(N), length(max_sizes));
+total_times = fill(NaN, length(N), length(max_sizes) + 1);
 ```
 
 ## Julia Timings
@@ -129,7 +129,6 @@ function time_open_modelica(n::Int)
 end
 
 function run_and_time_om!(ss_times, times, max_sizes, i, n)
-  run_and_time_julia!(ss_times, times, max_sizes, i, n)
   if n <= max_sizes[8]
     total_times[i, end] = time_open_modelica(n)
   end
@@ -197,11 +196,6 @@ for (i, timecat) in enumerate(("ODEProblem + f!", "Run", "Solve"))
   end
   Legend(f[i+1, 2], _lines, method_names)
 end
-let method_names_m = vcat(method_names, "OpenModelica");
-  ax = Axis(f[5, 1]; yscale = log10, xscale = log10, title = "Total Time")
-  _lines = map(Base.Fix1(lines!, N), eachcol(total_times))
-  Legend(f[5, 2], _lines, method_names_m)
-end
 f
 ```
 
@@ -218,7 +212,6 @@ Legend(f2[1,2], _lines, names)
 f2
 ```
 
-
 All three backends compiled more quickly with loops, but the C and LLVM backends are so much quicker to compile than the Julia backend that this made much less difference for them.
 The impact on runtime was more varied.
 
diff --git a/benchmarks/ModelingToolkit/ThermalFluid.jmd b/benchmarks/ModelingToolkit/ThermalFluid.jmd
index 2a03df05b..f46547ef8 100644
--- a/benchmarks/ModelingToolkit/ThermalFluid.jmd
+++ b/benchmarks/ModelingToolkit/ThermalFluid.jmd
@@ -302,12 +302,15 @@ end
 ```julia
 N = [5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800, 960, 1280];
 N_states = 4 .* N; # x-axis for plots
-# max size we test per method
-max_sizes = [480, last(N), last(N), last(N), last(N), last(N)];
+# max size we run per method, including Julia methods and OpenModelica
+# we do not run Dymola in this script, therefore it is excluded
+max_sizes = [480, last(N), last(N), last(N), last(N)];
 # NaN-initialize so Makie will ignore incomplete
 ss_times = fill(NaN, length(N), 2);
+# more indepth Julia-times
 times = fill((NaN,NaN,NaN), length(N), length(max_sizes) - 1);
-total_times = fill(NaN, length(N), length(max_sizes)+1); # +1 for Dymola
+# all times, including Dymola
+total_times = fill(NaN, length(N), length(max_sizes) + 1); # +1 for Dymola
 ```
 
 ## Julia Timings
@@ -374,7 +377,7 @@ total_times[:, 6] = translation_and_total_times[1:length(N),2]
 
 ```julia
 f = Figure(size=(800,1200));
-ss_names = ["MTK", "JSIR-Scalar", "JSIR-Loop"];
+ss_names = ["MTK", "JSIR"];
 let ax = Axis(f[1, 1]; yscale = log10, xscale = log10, title="Structural Simplify Time")
   _lines = map(eachcol(ss_times)) do ts
     lines!(N, ts)