diff --git a/src/IntervalOptimisation.jl b/src/IntervalOptimisation.jl
index 65cd52c..3de776a 100644
--- a/src/IntervalOptimisation.jl
+++ b/src/IntervalOptimisation.jl
@@ -2,16 +2,21 @@
 
 module IntervalOptimisation
 
+using IntervalArithmetic, IntervalRootFinding
+using LinearAlgebra
+
 export minimise, maximise,
        minimize, maximize
 
+export mean_value_form_scalar, third_order_taylor_form_scalar
+
+
 include("SortedVectors.jl")
 using .SortedVectors
 
-using IntervalArithmetic, IntervalRootFinding
-
 
 include("optimise.jl")
+include("centered_forms.jl")
 
 const minimize = minimise
 const maximize = maximise
diff --git a/src/centered_forms.jl b/src/centered_forms.jl
new file mode 100644
index 0000000..8116503
--- /dev/null
+++ b/src/centered_forms.jl
@@ -0,0 +1,34 @@
+
+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)
+
+"""
+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.v) ⋅ (X - m)
+end
+
+mean_value_form_scalar(f) = X -> mean_value_form_scalar(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 = 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_scalar(f) = X -> third_order_taylor_form_scalar(f, X)
diff --git a/src/optimise.jl b/src/optimise.jl
index 4629811..659b4c3 100644
--- a/src/optimise.jl
+++ b/src/optimise.jl
@@ -1,4 +1,7 @@
 
+interval_mid(X::Interval) = Interval(mid(X))
+interval_mid(X::IntervalBox) = IntervalBox(mid(X))
+
 """
     minimise(f, X, tol=1e-3)
 
@@ -27,7 +30,7 @@ function minimise(f, X::T, tol=1e-3) where {T}
         end
 
         # find candidate for upper bound of global minimum by just evaluating a point in the interval:
-        m = sup(f(Interval.(mid.(X))))   # evaluate at midpoint of current interval
+        m = sup(f(interval_mid(X)))   # evaluate at midpoint of current interval
 
         if m < global_min
             global_min = m