diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 80b0a51f..f5d11969 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -24,6 +24,7 @@ jobs:
     strategy:
       matrix:
         version:
+          - 'pre'
           - '1.10'
           # - '1.9'
           # - 'nightly'
diff --git a/Project.toml b/Project.toml
index 100adb69..ef7f5e72 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>", "Orjan Ameye <orjan.ameye@hotmail.com>"]
-version = "0.10.0"
+version = "0.10.1"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
@@ -17,6 +17,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Peaks = "18e31ff7-3703-566c-8e60-38913d67486b"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -39,13 +40,13 @@ BijectiveHilbert = "0.3.0"
 DSP = "0.7.9"
 DelimitedFiles = "1.9"
 Distances = "0.10.11"
-Documenter = "1.4"
 DocStringExtensions = "0.9.3"
+Documenter = "1.4"
 ExplicitImports = "1.6"
 FFTW = "1.8"
 HomotopyContinuation = "2.9"
-JLD2 = "0.4.48"
 JET = "0.9"
+JLD2 = "0.4.48"
 Latexify = "0.16"
 ModelingToolkit = "9.17"
 NonlinearSolve = "3.12"
@@ -53,18 +54,18 @@ OrderedCollections = "1.6"
 OrdinaryDiffEq = "v6.82"
 Peaks = "0.5"
 Plots = "1.39"
+PrecompileTools = "1.2"
 ProgressMeter = "1.7.2"
 SteadyStateDiffEq = "2"
 SymbolicUtils = "2.0"
 Symbolics = "5.30"
-
 julia = "1.10.0"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
diff --git a/benchmark/parametron.jl b/benchmark/parametron.jl
new file mode 100644
index 00000000..4b257a8c
--- /dev/null
+++ b/benchmark/parametron.jl
@@ -0,0 +1,25 @@
+using HarmonicBalance
+
+using Random
+const SEED = 0xd8e5d8df
+Random.seed!(SEED)
+
+@variables Ω γ λ F x θ η α ω0 ω t T ψ
+@variables x(t)
+
+natural_equation =
+    d(d(x, t), t) +
+    γ * d(x, t) +
+    Ω^2 * (1 - λ * cos(2 * ω * t + ψ)) * x +
+    α * x^3 +
+    η * d(x, t) * x^2
+forces = F * cos(ω * t + θ)
+dEOM = DifferentialEquation(natural_equation + forces, x)
+add_harmonic!(dEOM, x, ω)
+
+@time harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t);
+@time prob = HarmonicBalance.Problem(harmonic_eq)
+
+fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0)
+varied = ω => range(0.9, 1.1, 20)
+@time res = get_steady_states(prob, varied, fixed; show_progress=false)
diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl
index fbb47088..09fb529f 100644
--- a/src/HarmonicBalance.jl
+++ b/src/HarmonicBalance.jl
@@ -16,6 +16,8 @@ const HC = HomotopyContinuation
 using Plots: Plots, plot, plot!, savefig, heatmap, Plot
 using Latexify: Latexify, latexify
 
+using PrecompileTools: @setup_workload, @compile_workload
+
 # default global settings
 IM_TOL::Float64 = 1E-6
 function set_imaginary_tolerance(x::Float64)
@@ -71,4 +73,14 @@ export get_krylov_equations
 include("modules/FFTWExt.jl")
 using .FFTWExt
 
+@setup_workload begin
+    # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
+    # precompile file and potentially make loading faster.
+    @compile_workload begin
+        # all calls in this block will be precompiled, regardless of whether
+        # they belong to your package or not (on Julia 1.8 and higher)
+        include("precompilation.jl")
+    end
+end
+
 end # module
diff --git a/src/precompilation.jl b/src/precompilation.jl
new file mode 100644
index 00000000..dd269803
--- /dev/null
+++ b/src/precompilation.jl
@@ -0,0 +1,19 @@
+@variables Ω γ λ F x θ η α ω0 ω t T ψ
+@variables x(t)
+
+natural_equation =
+    d(d(x, t), t) +
+    γ * d(x, t) +
+    Ω^2 * (1 - λ * cos(2 * ω * t + ψ)) * x +
+    α * x^3 +
+    η * d(x, t) * x^2
+forces = F * cos(ω * t + θ)
+dEOM = DifferentialEquation(natural_equation + forces, x)
+add_harmonic!(dEOM, x, ω)
+
+harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t);
+prob = HarmonicBalance.Problem(harmonic_eq)
+
+fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0)
+varied = ω => range(0.9, 1.1, 20)
+res = get_steady_states(prob, varied, fixed; show_progress=false)