diff --git a/Project.toml b/Project.toml
index 156fa26..ef83540 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,13 +5,15 @@ version = "1.2.0"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-julia = "1"
 CommonSolve = "0.2"
+PrecompileTools = "1"
 SciMLBase = "2"
+julia = "1"
 
 [extras]
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
diff --git a/src/FiniteVolumeMethod1D.jl b/src/FiniteVolumeMethod1D.jl
index f8705c6..04fa887 100644
--- a/src/FiniteVolumeMethod1D.jl
+++ b/src/FiniteVolumeMethod1D.jl
@@ -25,4 +25,6 @@ include("solve.jl")
     (obj::Returns)(@nospecialize(args...); @nospecialize(kw...)) = obj.value
 end
 
+include("precompilation.jl")
+
 end # module
diff --git a/src/precompilation.jl b/src/precompilation.jl
new file mode 100644
index 0000000..1179e27
--- /dev/null
+++ b/src/precompilation.jl
@@ -0,0 +1,39 @@
+using PrecompileTools
+
+@setup_workload begin
+    # Minimal setup - these are lightweight and don't add much to precompilation time
+    mesh_points = LinRange(0.0, 1.0, 10)
+
+    @compile_workload begin
+        # Precompile FVMGeometry creation - common entry point
+        geom = FVMGeometry(mesh_points)
+
+        # Precompile boundary condition types with common cases
+        lhs_dir = Dirichlet(0.0)
+        rhs_dir = Dirichlet(1.0)
+        lhs_neu = Neumann(0.0)
+        rhs_neu = Neumann(0.0)
+
+        # Precompile BoundaryConditions with different combinations
+        bc_dir = BoundaryConditions(lhs = lhs_dir, rhs = rhs_dir)
+        bc_neu = BoundaryConditions(lhs = lhs_neu, rhs = rhs_neu)
+        bc_mix = BoundaryConditions(lhs = lhs_dir, rhs = rhs_neu)
+
+        # Precompile FVMProblem creation
+        diffusion_function = (u, x, t, p) -> one(u)
+        initial_condition = collect(mesh_points)
+
+        prob = FVMProblem(
+            mesh_points, lhs_dir, rhs_dir;
+            diffusion_function = diffusion_function,
+            initial_condition = initial_condition,
+            final_time = 1.0
+        )
+
+        # Precompile ODEProblem creation - this is a common operation
+        ode_prob = SciMLBase.ODEProblem(prob)
+
+        # Precompile jacobian_sparsity
+        jac_sp = jacobian_sparsity(prob)
+    end
+end