diff --git a/Project.toml b/Project.toml
index 4847e9e..0e8660b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "HarmonicSteadyState"
 uuid = "1158f75c-a779-4b85-8bfb-8fcf6bf02ced"
-version = "0.5.1"
+version = "0.5.2"
 authors = ["Orjan Ameye <orjan.ameye@hotmail.com>", "Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
 
 [deps]
@@ -75,7 +75,7 @@ SymbolicUtils = "4"
 Symbolics = "7"
 Test = "1.10"
 TestExtras = "0.3"
-julia = "1.12"
+julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
@@ -100,8 +100,5 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 
-[sources]
-HarmonicBalance = {url = "https://github.com/QuantumEngineeredSystems/HarmonicBalance.jl", rev = "source"}
-
 [targets]
 test = ["Aqua", "JET", "CheckConcreteStructs", "SymbolicUtils", "Symbolics", "ExplicitImports", "Test", "Random", "OrdinaryDiffEqTsit5", "Documenter", "Plots", "SteadyStateDiffEq", "LinearSolve", "TestExtras", "ModelingToolkitBase", "OrdinaryDiffEqRosenbrock", "NonlinearSolve", "HarmonicBalance", "Peaks", "Latexify", "QuantumCumulants"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 1d17185..006fa0b 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -15,4 +15,4 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 [compat]
 Documenter = "1"
 QuestBase = "0.4.1"
-julia = "1.12"
+julia = "1.10"
diff --git a/src/HarmonicSteadyState.jl b/src/HarmonicSteadyState.jl
index bcdc97b..0ff560d 100644
--- a/src/HarmonicSteadyState.jl
+++ b/src/HarmonicSteadyState.jl
@@ -27,8 +27,8 @@ using Random: Random # for setting seed
 using Distances: Distances
 using BijectiveHilbert: BijectiveHilbert, Simple2D, decode_hilbert!, encode_hilbert
 using HomotopyContinuation: HomotopyContinuation
-using Symbolics: Symbolics, unwrap, wrap, Num, get_variables
-using SymbolicUtils: SymbolicUtils
+using Symbolics: Symbolics, wrap, Num, get_variables
+using SymbolicUtils: SymbolicUtils, unwrap
 
 const HC = HomotopyContinuation
 import FunctionWrappers: FunctionWrapper
diff --git a/src/LinearResponse/Lorentzian_spectrum.jl b/src/LinearResponse/Lorentzian_spectrum.jl
index 570bfcc..ae939be 100644
--- a/src/LinearResponse/Lorentzian_spectrum.jl
+++ b/src/LinearResponse/Lorentzian_spectrum.jl
@@ -114,7 +114,7 @@ function JacobianSpectrum(
     uv_ωnums = [
         real(
             SymbolicUtils.unwrap_const(
-                Symbolics.unwrap(Symbolics.substitute(hvars[pair][1].ω, substitutions))
+                SymbolicUtils.unwrap(Symbolics.substitute(hvars[pair][1].ω, substitutions))
             ),
         ) for pair in uv_pairs
     ]
diff --git a/src/Result.jl b/src/Result.jl
index 3ec3a10..9a883b2 100644
--- a/src/Result.jl
+++ b/src/Result.jl
@@ -65,7 +65,7 @@ function Result(
     )
 end
 
-Symbolics.get_variables(res::Result)::Vector{Num} = get_variables(res.problem)
+Symbolics.get_variables(res::Result) = get_variables(res.problem)
 
 function Base.show(io::IO, r::Result)
     println(io, "A steady state result for ", length(r.solutions), " parameter points")
diff --git a/src/transform_solutions.jl b/src/transform_solutions.jl
index 17409cc..568f4bd 100644
--- a/src/transform_solutions.jl
+++ b/src/transform_solutions.jl
@@ -310,12 +310,14 @@ function _to_lab_frame_velocity(soln, vars, times)
     for var in vars
         val = real(
             SymbolicUtils.unwrap_const(
-                Symbolics.unwrap(
-                    substitute_all(Symbolics.unwrap(_remove_brackets(var)), soln)
+                SymbolicUtils.unwrap(
+                    substitute_all(SymbolicUtils.unwrap(_remove_brackets(var)), soln)
                 ),
             ),
         )
-        ω = real(SymbolicUtils.unwrap_const(Symbolics.unwrap(substitute_all(var.ω, soln))))
+        ω = real(
+            SymbolicUtils.unwrap_const(SymbolicUtils.unwrap(substitute_all(var.ω, soln)))
+        )
         if var.type == "u"
             timetrace .+= -ω * val * sin.(ω * times)
         elseif var.type == "v"
@@ -331,12 +333,14 @@ function _to_lab_frame(soln, vars, times)::Vector{AbstractFloat}
     for var in vars
         val = real(
             SymbolicUtils.unwrap_const(
-                Symbolics.unwrap(
-                    substitute_all(Symbolics.unwrap(_remove_brackets(var)), soln)
+                SymbolicUtils.unwrap(
+                    substitute_all(SymbolicUtils.unwrap(_remove_brackets(var)), soln)
                 ),
             ),
         )
-        ω = real(SymbolicUtils.unwrap_const(Symbolics.unwrap(substitute_all(var.ω, soln))))
+        ω = real(
+            SymbolicUtils.unwrap_const(SymbolicUtils.unwrap(substitute_all(var.ω, soln)))
+        )
         if var.type == "u"
             timetrace .+= val * cos.(ω * times)
         elseif var.type == "v"