diff --git a/Project.toml b/Project.toml
index 523c7b3..2ff7d09 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,6 +13,7 @@ FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
@@ -32,6 +33,7 @@ FastTransforms = "0.12, 0.13, 0.14, 0.15.1"
 FillArrays = "0.11, 0.12, 0.13, 1"
 IntervalSets = "0.5, 0.6, 0.7"
 LazyArrays = "0.22, 1"
+PrecompileTools = "1"
 Reexport = "0.2, 1"
 SpecialFunctions = "0.10, 1.0, 2"
 Static = "0.8"
diff --git a/src/ApproxFunOrthogonalPolynomials.jl b/src/ApproxFunOrthogonalPolynomials.jl
index 9cd5ed5..a7907b7 100644
--- a/src/ApproxFunOrthogonalPolynomials.jl
+++ b/src/ApproxFunOrthogonalPolynomials.jl
@@ -142,5 +142,6 @@ include("specialfunctions.jl")
 include("fastops.jl")
 include("symeigen.jl")
 include("show.jl")
+include("precompile.jl")
 
 end
diff --git a/src/precompile.jl b/src/precompile.jl
new file mode 100644
index 0000000..0821764
--- /dev/null
+++ b/src/precompile.jl
@@ -0,0 +1,41 @@
+using PrecompileTools
+
+@setup_workload begin
+    splist = Any[Jacobi(1, 1), Chebyshev(), Ultraspherical(1)]
+    append!(splist, Any[Jacobi(1, 1, 0..1), Chebyshev(0..1), Ultraspherical(1, 0..1)])
+    spreal = copy(splist)
+    push!(splist, Chebyshev(Segment(1.0+im,2.0+2im)))
+    # special functions
+    spfns = (sin, cos, exp)
+    v = ones(2)
+    m = ones(2,2)
+    a = ones(2,2,2)
+    @compile_workload begin
+        for S in splist
+        	f = Fun(S,v)
+            f(0.1)
+            1/(f^2+1)
+        	abs(DefiniteIntegral()*f - sum(f))
+        	norm(Derivative()*f-f')
+        	norm(first(cumsum(f)))
+            for spfn in spfns
+                spfn(f)
+            end
+        end
+        for S in spreal
+            transform(S, v)
+            itransform(S, v)
+            transform!(S, v)
+            itransform!(S, v)
+            transform(S * S, m)
+            itransform(S * S, m)
+
+            [Derivative(S)^2 + 1; Dirichlet(S)] \ [1,0]
+
+            f2 = Fun((x,y) -> x*y, S^2); f2(0.1, 0.1)
+            P = ProductFun((x,y) -> x*y, S^2); P(0.1, 0.1)
+            L = LowRankFun((x,y) -> x*y, S^2); L(0.1, 0.1)
+        end
+        roots(Fun(sin,Interval(big"0.0", big"1.0")))
+    end
+end