diff --git a/Project.toml b/Project.toml
index 2aec45a..cd8f023 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.1.3"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -16,6 +17,7 @@ TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 DocStringExtensions = "0.8, 0.9"
 LinearAlgebra = "1"
 ModelingToolkit = "10.10"
+PrecompileTools = "1"
 RandomizedLinAlg = "0.1"
 Setfield = "0.8, 1"
 SparseArrays = "1"
diff --git a/src/ModelOrderReduction.jl b/src/ModelOrderReduction.jl
index 9691054..05be57d 100644
--- a/src/ModelOrderReduction.jl
+++ b/src/ModelOrderReduction.jl
@@ -19,4 +19,6 @@ export POD, reduce!
 include("deim.jl")
 export deim
 
+include("precompile.jl")
+
 end
diff --git a/src/precompile.jl b/src/precompile.jl
new file mode 100644
index 0000000..6777c0c
--- /dev/null
+++ b/src/precompile.jl
@@ -0,0 +1,35 @@
+using PrecompileTools
+
+@setup_workload begin
+    # Setup code - create minimal test data
+    n = 20  # small state dimension
+    m = 10  # small number of snapshots
+    snapshot_matrix = Float64[sin(i * j / n) for i in 1:n, j in 1:m]
+    snapshot_vov = [Float64[sin(i * j / n) for i in 1:n] for j in 1:m]
+
+    # Create an orthonormal basis for deim_interpolation_indices
+    Q, _ = qr(snapshot_matrix)
+    deim_basis = Matrix(Q[:, 1:5])
+
+    @compile_workload begin
+        # POD construction with matrix input
+        pod_mat = POD(snapshot_matrix, 3)
+
+        # POD construction with vector of vectors
+        pod_vov = POD(snapshot_vov, 3)
+
+        # reduce! with SVD algorithm
+        reduce!(pod_mat, SVD())
+
+        # reduce! with TSVD algorithm
+        pod_tsvd = POD(snapshot_matrix, 3)
+        reduce!(pod_tsvd, TSVD())
+
+        # reduce! with RSVD algorithm
+        pod_rsvd = POD(snapshot_matrix, 3)
+        reduce!(pod_rsvd, RSVD())
+
+        # deim_interpolation_indices (core DEIM algorithm)
+        deim_interpolation_indices(deim_basis)
+    end
+end