diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl
index d70ed783..4f5448dc 100644
--- a/src/IntervalRootFinding.jl
+++ b/src/IntervalRootFinding.jl
@@ -8,7 +8,7 @@ using IntervalArithmetic
 using ForwardDiff
 using StaticArrays
 
-using LinearAlgebra: I, Diagonal
+using LinearAlgebra: I, Diagonal, ⋅
 
 
 import Base: ⊆, show, big, \
@@ -24,7 +24,8 @@ export
     gauss_seidel_interval, gauss_seidel_interval!,
     gauss_seidel_contractor, gauss_seidel_contractor!,
     gauss_elimination_interval, gauss_elimination_interval!,
-    slope
+    slope,
+    mean_value_form_scalar, mean_value_form_vector, third_order_taylor_form_scalar
 
 export isunique, root_status
 
@@ -51,6 +52,7 @@ include("newton1d.jl")
 include("quadratic.jl")
 include("linear_eq.jl")
 include("slopes.jl")
+include("centered_forms.jl")
 
 
 gradient(f) = X -> ForwardDiff.gradient(f, X[:])
diff --git a/src/centered_forms.jl b/src/centered_forms.jl
new file mode 100644
index 00000000..90fe09f1
--- /dev/null
+++ b/src/centered_forms.jl
@@ -0,0 +1,54 @@
+import ForwardDiff: gradient, jacobian, hessian
+
+gradient(f, X::IntervalBox) = gradient(f, X.v)
+# jacobian(f, X::IntervalBox) = jacobian(f, X.v)
+hessian(f, X::IntervalBox) = hessian(f, X.v)
+
+import Base: *
+
+*(M::AbstractMatrix, X::IntervalBox) = IntervalBox(M * X.v)
+
+
+
+"""
+Calculate the mean-value form of a function ``f:\\mathbb{R}^n \\to \\mathbb{R}``
+using the gradient ``\nabla f``;
+this gives a second-order approximation.
+"""
+function mean_value_form_scalar(f, X)
+    m = IntervalBox(mid(X))
+
+    return f(m) + gradient(f, X) ⋅ (X - m)
+end
+
+mean_value_form_scalar(f) = X -> mean_value_form_scalar(f, X)
+
+"""
+Calculate the mean-value form of a function ``f:\\mathbb{R}^n \\to \\mathbb{R^n}``
+using the Jacobian;
+this gives a second-order approximation.
+"""
+function mean_value_form_vector(f, X)
+    m = IntervalBox(mid(X))
+
+    return f(m) + jacobian(f, X) * (X - m)
+end
+
+mean_value_form_vector(f) = X -> mean_value_form_vector(f, X)
+
+
+
+
+"""
+Calculate a third-order Taylor form of ``f:\\mathbb{R}^n \\to \\mathbb{R}`` using the Hessian.
+"""
+function third_order_taylor_form_scalar(f, X)
+    m = IntervalBox(mid(X))
+
+    H = ForwardDiff.hessian(f, X)
+    δ = X - m
+
+    return f(m) + gradient(f, m) ⋅ δ + 0.5 * sum(δ[i]*H[i,j]*δ[j] for i in 1:length(X) for j in 1:length(X))
+end
+
+third_order_taylor_form(f) = X -> third_order_taylor_form(f, X)
diff --git a/test/centered_forms.jl b/test/centered_forms.jl
new file mode 100644
index 00000000..2c8a00a0
--- /dev/null
+++ b/test/centered_forms.jl
@@ -0,0 +1,13 @@
+using IntervalArithmetic, IntervalRootFinding, StaticArrays
+using Test
+
+
+
+f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
+f(X) = f(X[1], X[2])
+X = (-6..6) × (-6..6)
+
+rts = roots(f, X, Bisection)
+rts2 = roots(mean_value_form_vector(f), X, Bisection)
+
+@test length(rts) == length(rts2) == 4
diff --git a/test/runtests.jl b/test/runtests.jl
index 148de3a7..eeda400b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,3 +9,4 @@ include("newton1d.jl")
 include("quadratic.jl")
 include("linear_eq.jl")
 include("slopes.jl")
+include("centered_forms.jl")