diff --git a/src/Symbolics.jl b/src/Symbolics.jl
index 2b6a45912..561311878 100644
--- a/src/Symbolics.jl
+++ b/src/Symbolics.jl
@@ -140,6 +140,9 @@ export symbolic_linear_solve, solve_for
 include("groebner_basis.jl")
 export groebner_basis, is_groebner_basis
 
+include("taylor.jl")
+export taylor
+
 import Libdl
 include("build_function.jl")
 export build_function
diff --git a/src/taylor.jl b/src/taylor.jl
new file mode 100644
index 000000000..c0fdf16f0
--- /dev/null
+++ b/src/taylor.jl
@@ -0,0 +1,31 @@
+# TODO: around arbitrary point x = x0 ≠ 0
+# TODO: error if x is not a "pure variable"
+# TODO: optimize for multiple orders with loop/recursion
+# TODO: get rational coefficients, not floats
+"""
+    taylor(f, x, n)
+
+Calculate the `n`-th order term(s) in the Taylor series of the expression `f(x)` around `x = 0`.
+
+Examples
+========
+```julia
+julia> @variables x
+1-element Vector{Num}:
+ x
+
+julia> taylor(exp(x), x, 0:3)
+1.0 + x + 0.5(x^2) + 0.16666666666666666(x^3)
+```
+"""
+function taylor(f, x, n::Int)
+    D = Differential(x)
+    n! = factorial(n)
+    c = (D^n)(f) / n!
+    c = expand_derivatives(c)
+    c = substitute(c, x => 0)
+    return c * x^n
+end
+function taylor(f, x, n::AbstractArray{Int})
+    return sum(taylor.(f, x, n))
+end
diff --git a/test/taylor.jl b/test/taylor.jl
new file mode 100644
index 000000000..4ce3edbe0
--- /dev/null
+++ b/test/taylor.jl
@@ -0,0 +1,25 @@
+@testset "Taylor series" begin
+    # https://en.wikipedia.org/wiki/Taylor_series#List_of_Maclaurin_series_of_some_common_functions
+    @variables x
+    @test taylor(exp(x), x, 0:9) - sum(1 / factorial(n) * x^n for n in 0:9) == 0
+    @test taylor(log(1-x), x, 0:9) - sum(-1/n * x^n for n in 1:9) == 0
+    @test taylor(log(1+x), x, 0:9) - sum((-1)^(n+1) * 1/n * x^n for n in 1:9) == 0
+
+    @test taylor(1/(1-x), x, 0:9) - sum(x^n for n in 0:9) == 0
+    @test taylor(1/(1-x)^2, x, 0:8) - sum(n * x^(n-1) for n in 1:9) == 0
+    @test taylor(1/(1-x)^3, x, 0:7) - sum((n-1)*n/2 * x^(n-2) for n in 2:9) == 0
+    for α in (-1//2, 0, 1//2, 1, 2, 3)
+        @test taylor((1+x)^α, x, 0:7) - sum(binomial(α, n) * x^n for n in 0:7) == 0
+    end
+
+    @test taylor(sin(x), x, 0:7) - sum((-1)^n/factorial(2*n+1) * x^(2*n+1) for n in 0:3) == 0
+    @test taylor(cos(x), x, 0:7) - sum((-1)^n/factorial(2*n) * x^(2*n) for n in 0:3) == 0
+    @test taylor(tan(x), x, 0:7) - taylor(taylor(sin(x), x, 0:7) / taylor(cos(x), x, 0:7), x, 0:7) == 0
+    @test taylor(asin(x), x, 0:7) - sum(factorial(2*n)/(4^n*factorial(n)^2*(2*n+1)) * x^(2*n+1) for n in 0:3) == 0
+    @test taylor(acos(x), x, 0:7) - taylor(π/2 - asin(x), x, 0:7) == 0
+    @test taylor(atan(x), x, 0:7) - taylor(asin(x/√(1+x^2)), x, 0:7) == 0
+
+    @test taylor(sinh(x), x, 0:7) - sum(1/factorial(2*n+1) * x^(2*n+1) for n in 0:3) == 0
+    @test taylor(cosh(x), x, 0:7) - sum(1/factorial(2*n) * x^(2*n) for n in 0:3) == 0
+    @test taylor(tanh(x), x, 0:7) - (x - x^3/3 + 2/15*x^5 - 17/315*x^7) == 0
+end