diff --git a/Project.toml b/Project.toml
index d6c26eb..ae8e87b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Richardson = "708f8203-808e-40c0-ba2d-98a6953ed40d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 ChainRulesCore = "0.10.3, 1"
@@ -19,6 +20,7 @@ Random = "<0.0.1, 1"
 Richardson = "1.2"
 SparseArrays = "<0.0.1, 1"
 StaticArrays = "0.12, 1.0"
+Unitful = "1.24.0"
 julia = "1"
 
 [extras]
diff --git a/src/FiniteDifferences.jl b/src/FiniteDifferences.jl
index 82dd64c..21f8cd3 100644
--- a/src/FiniteDifferences.jl
+++ b/src/FiniteDifferences.jl
@@ -7,6 +7,7 @@ using Random
 using Richardson
 using SparseArrays
 using StaticArrays
+using Unitful
 
 export to_vec, grad, jacobian, jvp, j′vp
 
diff --git a/src/methods.jl b/src/methods.jl
index f3d8ea7..63d6314 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -107,9 +107,9 @@ end
     FiniteDifferenceMethod(
         grid::AbstractVector{Int},
         q::Int;
-        condition::Real=DEFAULT_CONDITION,
-        factor::Real=DEFAULT_FACTOR,
-        max_range::Real=Inf
+        condition::Number=DEFAULT_CONDITION,
+        factor::Number=DEFAULT_FACTOR,
+        max_range::Number=Inf
     )
 
 Construct a finite difference method.
@@ -120,9 +120,9 @@ Construct a finite difference method.
 - `q::Int`: Order of the derivative to estimate.
 
 # Keywords
-- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref).
-- `factor::Real`: Factor number. See [`DEFAULT_FACTOR`](@ref).
-- `max_range::Real=Inf`: Maximum distance that a function is evaluated from the input at
+- `condition::Number`: Condition number. See [`DEFAULT_CONDITION`](@ref).
+- `factor::Number`: Factor number. See [`DEFAULT_FACTOR`](@ref).
+- `max_range::Number=Inf`: Maximum distance that a function is evaluated from the input at
     which the derivative is estimated.
 
 # Returns
@@ -131,9 +131,9 @@ Construct a finite difference method.
 function FiniteDifferenceMethod(
     grid::SVector{P,Int},
     q::Int;
-    condition::Real=DEFAULT_CONDITION,
-    factor::Real=DEFAULT_FACTOR,
-    max_range::Real=Inf,
+    condition::Number=DEFAULT_CONDITION,
+    factor::Number=DEFAULT_FACTOR,
+    max_range::Number=Inf,
 ) where P
     _check_p_q(P, q)
     coefs, coefs_neighbourhood, ∇f_magnitude_mult, f_error_mult = _coefs_mults(grid, q)
@@ -153,7 +153,7 @@ function FiniteDifferenceMethod(grid::AbstractVector{Int}, q::Int; kw_args...)
 end
 
 """
-    (m::FiniteDifferenceMethod)(f, x::T) where T<:AbstractFloat
+    (m::FiniteDifferenceMethod)(f, x::T) where T<:Number
 
 Estimate the derivative of `f` at `x` using the finite differencing method `m` and an
 automatically determined step size.
@@ -188,12 +188,12 @@ julia> FiniteDifferences.estimate_step(fdm, sin, 1.0)  # Computes step size and
 # We loop over all concrete subtypes of `FiniteDifferenceMethod` for Julia v1.0 compatibility.
 for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
     @eval begin
-        function (m::$T)(f::TF, x::Real) where TF
+        function (m::$T)(f::TF, x::Number) where TF
             x = float(x)  # Assume that converting to float is desired, if it isn't already.
             step = first(estimate_step(m, f, x))
             return m(f, x, step)
         end
-        function (m::$T{P,0})(f::TF, x::Real) where {P,TF}
+        function (m::$T{P,0})(f::TF, x::Number) where {P,TF}
             # The automatic step size calculation fails if `Q == 0`, so handle that edge
             # case.
             return f(x)
@@ -202,7 +202,7 @@ for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
 end
 
 """
-    (m::FiniteDifferenceMethod)(f, x::T, step::Real) where T<:AbstractFloat
+    (m::FiniteDifferenceMethod)(f, x::T, step::Number) where T<:Number
 
 Estimate the derivative of `f` at `x` using the finite differencing method `m` and a given
 step size.
@@ -210,7 +210,7 @@ step size.
 # Arguments
 - `f`: Function to estimate derivative of.
 - `x::T`: Input to estimate derivative at.
-- `step::Real`: Step size.
+- `step::Number`: Step size.
 
 # Returns
 - Estimate of the derivative.
@@ -235,7 +235,7 @@ julia> fdm(sin, 1, 1e-3) - cos(1)  # Check the error.
 # We loop over all concrete subtypes of `FiniteDifferenceMethod` for 1.0 compatibility.
 for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
     @eval begin
-        function (m::$T{P,Q})(f::TF, x::Real, step::Real) where {P,Q,TF}
+        function (m::$T{P,Q})(f::TF, x::Number, step::Number) where {P,Q,TF}
             x = float(x)  # Assume that converting to float is desired, if it isn't already.
             fs = _eval_function(m, f, x, step)
             return _compute_estimate(m, fs, x, step, m.coefs)
@@ -244,8 +244,8 @@ for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
 end
 
 function _eval_function(
-    m::FiniteDifferenceMethod, f::TF, x::T, step::Real,
-) where {TF,T<:AbstractFloat}
+    m::FiniteDifferenceMethod, f::TF, x::T, step::Number,
+) where {TF,T<:Number}
     return f.(x .+ T(step) .* m.grid)
 end
 
@@ -253,14 +253,17 @@ function _compute_estimate(
     m::FiniteDifferenceMethod{P,Q},
     fs::SVector{P,TF},
     x::T,
-    step::Real,
+    step::Number,
     coefs::SVector{P,Float64},
-) where {P,Q,TF,T<:AbstractFloat}
+) where {P,Q,TF,T<:Number}
     # If we substitute `T.(coefs)` in the expression below, then allocations occur. We
     # therefore perform the broadcasting first. See
     # https://github.com/JuliaLang/julia/issues/39151.
-    _coefs = T.(coefs)
-    return sum(fs .* _coefs) ./ T(step)^Q
+    #
+    # We strip units because the estimate coefficients are just weights for values of f.
+    N = numType(T)
+    _coefs = N.(coefs)
+    return sum(fs .* _coefs) ./ T(step) ^ Q
 end
 
 # Check the method and derivative orders for consistency.
@@ -338,7 +341,7 @@ end
         m::FiniteDifferenceMethod,
         f,
         x::T
-    ) where T<:AbstractFloat
+    ) where T<:Number
 
 Estimate the step size for a finite difference method `m`. Also estimates the error of the
 estimate of the derivative.
@@ -349,31 +352,39 @@ estimate of the derivative.
 - `x::T`: Point to estimate the derivative at.
 
 # Returns
-- `Tuple{<:AbstractFloat, <:AbstractFloat}`: Estimated step size and an estimate of the
+- `Tuple{<:Number, <:Number}`: Estimated step size and an estimate of the
     error of the finite difference estimate. The error will be `NaN` if the method failed
     to estimate the error.
 """
 function estimate_step(
     m::UnadaptedFiniteDifferenceMethod, f::TF, x::T,
-) where {TF,T<:AbstractFloat}
+) where {TF,T<:Number}
     step, acc = _compute_step_acc_default(m, x)
+    xunit = unit(x)
+    step = withUnit(xunit, step)
+    acc = withUnit(xunit,acc)
     return _limit_step(m, x, step, acc)
 end
 function estimate_step(
     m::AdaptedFiniteDifferenceMethod{P,Q}, f::TF, x::T,
-) where {P,Q,TF,T<:AbstractFloat}
+) where {P,Q,TF,T<:Number}
     ∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
-    if ∇f_magnitude == 0.0 || f_magnitude == 0.0
-        step, acc = _compute_step_acc_default(m, x)
-    else
-        step, acc = _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))
-    end
+    xunit = unit(x)
+    dfunit = unit(first(f(x))) / unit(x) ^ Q
+    step, acc =
+        if ∇f_magnitude == withUnit(unit(∇f_magnitude),0.0) || f_magnitude == withUnit(unit(f_magnitude), 0.0)
+            _compute_step_acc_default(m, x)
+        else
+            _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))
+        end
+    step = withUnit(xunit, step)
+    acc = withUnit(dfunit, acc)
     return _limit_step(m, x, step, acc)
 end
 
 function _estimate_magnitudes(
     m::FiniteDifferenceMethod{P,Q}, f::TF, x::T,
-) where {P,Q,TF,T<:AbstractFloat}
+) where {P,Q,TF,T<:Number}
     step = first(estimate_step(m, f, x))
     fs = _eval_function(m, f, x, step)
     # Estimate magnitude of `∇f` in a neighbourhood of `x`.
@@ -388,39 +399,44 @@ function _estimate_magnitudes(
     return ∇f_magnitude, f_magnitude
 end
 
-function _compute_step_acc_default(m::FiniteDifferenceMethod, x::T) where {T<:AbstractFloat}
+function _compute_step_acc_default(m::FiniteDifferenceMethod, x::T) where {T<:Number}
     # Compute a default step size using a heuristic and [`DEFAULT_CONDITION`](@ref).
     return _compute_step_acc(m, m.condition, eps(T))
 end
 
 function _compute_step_acc(
-    m::FiniteDifferenceMethod{P,Q}, ∇f_magnitude::Real, f_error::Real,
+    m::FiniteDifferenceMethod{P,Q}, ∇f_magnitude::Number, f_error::Number,
 ) where {P,Q}
     # Set the step size by minimising an upper bound on the error of the estimate.
     C₁ = f_error * m.f_error_mult * m.factor
     C₂ = ∇f_magnitude * m.∇f_magnitude_mult
-    step = (Q / (P - Q) * (C₁ / C₂))^(1 / P)
+    step = (Q / (P - Q) * (C₁ / C₂))^(1 // P)
     # Estimate the accuracy of the method.
     acc = C₁ * step^(-Q) + C₂ * step^(P - Q)
     return step, acc
 end
 
 function _limit_step(
-    m::FiniteDifferenceMethod, x::T, step::Real, acc::Real,
-) where {T<:AbstractFloat}
+    m::FiniteDifferenceMethod, x::T, step::Number, acc::Number,
+) where {T<:Number}
+    xunit = unit(x)
     # First, limit the step size based on the maximum range.
-    step_max = m.max_range / maximum(abs.(m.grid))
+    step_max = withUnit(
+        xunit,
+        m.max_range / maximum(abs.(m.grid))
+    )
     if step > step_max
         step = step_max
-        acc = NaN
+        acc = withUnit(xunit,NaN)
     end
     # Second, prevent very large step sizes, which can occur for high-order methods or
     # slowly-varying functions.
     step_default, _ = _compute_step_acc_default(m, x)
+    step_default = withUnit(xunit, step_default)
     step_max_default = 1000step_default
     if step > step_max_default
         step = step_max_default
-        acc = NaN
+        acc = withUnit(xunit,NaN)
     end
     return step, acc
 end
@@ -433,9 +449,9 @@ for direction in [:forward, :central, :backward]
             p::Int,
             q::Int;
             adapt::Int=1,
-            condition::Real=DEFAULT_CONDITION,
-            factor::Real=DEFAULT_FACTOR,
-            max_range::Real=Inf,
+            condition::Number=DEFAULT_CONDITION,
+            factor::Number=DEFAULT_FACTOR,
+            max_range::Number=Inf,
             geom::Bool=false
         )
             _check_p_q(p, q)
@@ -484,9 +500,9 @@ for direction in [:forward, :central, :backward]
         p::Int,
         q::Int;
         adapt::Int=1,
-        condition::Real=DEFAULT_CONDITION,
-        factor::Real=DEFAULT_FACTOR,
-        max_range::Real=Inf,
+        condition::Number=DEFAULT_CONDITION,
+        factor::Number=DEFAULT_FACTOR,
+        max_range::Number=Inf,
         geom::Bool=false
     )
 
@@ -500,9 +516,9 @@ Construct a finite difference method at a $($(Meta.quot(direction))) grid of `p`
 - `adapt::Int=1`: Use another finite difference method to estimate the magnitude of the
     `p`th order derivative, which is important for the step size computation. Recurse
     this procedure `adapt` times.
-- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref).
-- `factor::Real`: Factor number. See [`DEFAULT_FACTOR`](@ref).
-- `max_range::Real=Inf`: Maximum distance that a function is evaluated from the input at
+- `condition::Number`: Condition number. See [`DEFAULT_CONDITION`](@ref).
+- `factor::Number`: Factor number. See [`DEFAULT_FACTOR`](@ref).
+- `max_range::Number=Inf`: Maximum distance that a function is evaluated from the input at
     which the derivative is estimated.
 - `geom::Bool`: Use geometrically spaced points instead of linearly spaced points.
 
@@ -552,10 +568,10 @@ end
     extrapolate_fdm(
         m::FiniteDifferenceMethod,
         f,
-        x::Real,
-        initial_step::Real=10,
+        x::Number,
+        initial_step::Number=10,
         power::Int=1,
-        breaktol::Real=Inf,
+        breaktol::Number=Inf,
         kw_args...
     )
 
@@ -568,19 +584,19 @@ automatically sets `power = 2` if `m` is symmetric and `power = 1`. Moreover, it
 # Arguments
 - `m::FiniteDifferenceMethod`: Finite difference method to estimate the step size for.
 - `f`: Function to evaluate the derivative of.
-- `x::Real`: Point to estimate the derivative at.
-- `initial_step::Real=10`: Initial step size.
+- `x::Number`: Point to estimate the derivative at.
+- `initial_step::Number=10`: Initial step size.
 
 # Returns
-- `Tuple{<:AbstractFloat, <:AbstractFloat}`: Estimate of the derivative and error.
+- `Tuple{<:Number, <:Number}`: Estimate of the derivative and error.
 """
 function extrapolate_fdm(
     m::FiniteDifferenceMethod,
     f,
-    x::Real,
-    initial_step::Real=10;
+    x::Number,
+    initial_step::Number=10;
     power::Int=1,
-    breaktol::Real=Inf,
+    breaktol::Number=Inf,
     kw_args...
 )
     (power == 1 && _is_symmetric(m)) && (power = 2)
@@ -592,3 +608,40 @@ function extrapolate_fdm(
         kw_args...
     )
 end
+
+"""
+Attaches a given unit to a given value, if it is dimensionless.
+If the value is not dimensionless, attempts a conversion to the given unit.
+"""
+function withUnit(targetUnit, value)
+
+    if Unitful.dimension(value) == Unitful.NoDims
+
+        value * targetUnit
+
+    else
+
+        Unitful.uconvert(targetUnit, value)
+
+    end # if
+
+end # function
+
+"""
+Retrieves the number type of a quantity, or returns the type itself in the case of a raw number.
+"""
+function numType(x::Number)
+    typeof(x)
+end # function
+
+function numType(x::Type{<:Number})
+    x
+end
+
+function numType(x::Unitful.AbstractQuantity)
+    Unitful.numtype(typeof(x))
+end # function
+
+function numType(x::Type{<:Unitful.AbstractQuantity})
+    Unitful.numtype(x)
+end
diff --git a/test/methods.jl b/test/methods.jl
index af8e15b..37f6f4c 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -82,6 +82,8 @@ struct NotAFunction end # not <: Function on purpose, cf #224
     @testset "Test allocations" begin
         m = central_fdm(5, 2, adapt=2)
         @test @ballocated($m(sin, 1)) == 0
+        fn(t) = sin(t / u"s")
+        @test @ballocated($m($fn, 1u"s")) == 0
     end
 
     # Integration test to ensure that Integer-output functions can be tested.
@@ -89,6 +91,21 @@ struct NotAFunction end # not <: Function on purpose, cf #224
         @test isapprox(central_fdm(5, 1)(x -> 5, 0), 0; rtol=1e-12, atol=1e-12)
     end
 
+    # Integration test to ensure that Unitful-output functions can be tested.
+    @testset "Unitful output" begin
+
+        fn(x) = 5.0u"J/s" * x
+
+        derivativeVal = 5.0u"J/s"
+
+        adaptedDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s",1)
+
+        @test unit(adaptedDerivativeVal) == u"J/s"
+
+        @test isapprox(adaptedDerivativeVal, derivativeVal; rtol=1e-12, atol=1e-12u"J/s")
+
+    end
+
     @testset "Adaptation improves estimate" begin
         @test abs(forward_fdm(5, 1, adapt=0)(log, 0.001) - 1000) > 1
         @test forward_fdm(5, 1, adapt=1)(log, 0.001) ≈ 1000
diff --git a/test/runtests.jl b/test/runtests.jl
index 63039f7..a2d524f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,7 @@ using Random
 using SparseArrays
 using StaticArrays
 using Test
+using Unitful
 
 
 Random.seed!(1)