diff --git a/src/Interpolations.jl b/src/Interpolations.jl
index b7cd1e3e..91a50669 100644
--- a/src/Interpolations.jl
+++ b/src/Interpolations.jl
@@ -28,6 +28,7 @@ export
     # see the following files for further exports:
     # b-splines/b-splines.jl
     # extrapolation/extrapolation.jl
+    # monotonic/monotonic.jl
     # scaling/scaling.jl
 
 using LinearAlgebra, SparseArrays
@@ -388,6 +389,7 @@ end
 include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
 include("gridded/gridded.jl")
+include("monotonic/monotonic.jl")
 include("extrapolation/extrapolation.jl")
 include("scaling/scaling.jl")
 include("utils.jl")
diff --git a/src/io.jl b/src/io.jl
index da0bdffb..eae584f2 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -44,6 +44,20 @@ function Base.showarg(io::IO, A::ScaledInterpolation{T}, toplevel) where {T}
     end
 end
 
+function Base.showarg(io::IO, A::MonotonicInterpolation{T, TCoeffs, Tel, Type, K}, toplevel) where {T, TCoeffs, Tel, Type, K}
+    print(io, "interpolate(")
+    _showknots(io, A.knots)
+    print(io, ", ")
+    Base.showarg(io, A.A, false)
+    print(io, ", ")
+    show(io, A.it)
+    if toplevel
+        print(io, ") with element type ",T)
+    else
+        print(io, ')')
+    end
+end
+
 function Base.showarg(io::IO, A::Extrapolation{T,N,TI,IT,ET}, toplevel) where {T,N,TI,IT,ET}
     print(io, "extrapolate(")
     Base.showarg(io, A.itp, false)
diff --git a/src/monotonic/monotonic.jl b/src/monotonic/monotonic.jl
new file mode 100644
index 00000000..fb385851
--- /dev/null
+++ b/src/monotonic/monotonic.jl
@@ -0,0 +1,373 @@
+
+#=
+    Interpolate using cubic Hermite splines. The breakpoints in arrays xbp and ybp are assumed to be sorted.
+    Evaluate the function in all points of the array xeval.
+    Methods:
+        "Linear"                yuck
+        "FiniteDifference"      classic cubic interpolation, no tension parameter
+                                Finite difference can overshoot for non-monotonic data
+        "Cardinal"              cubic cardinal splines, uses tension parameter which must be between [0,1]
+                                cubin cardinal splines can overshoot for non-monotonic data
+                                (increasing tension decreases overshoot)
+        "FritschCarlson"        monotonic - tangents are first initialized, then adjusted if they are not monotonic
+                                can overshoot for non-monotonic data
+        "FritschButland"        monotonic - faster algorithm (only requires one pass) but somewhat higher apparent "tension"
+        "Steffen"               monotonic - also only one pass, results usually between FritschCarlson and FritschButland
+    Sources:
+        Fritsch & Carlson (1980), "Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.
+        Fritsch & Butland (1984), "A Method for Constructing Local Monotone Piecewise Cubic Interpolants", doi:10.1137/0905021.
+        Steffen (1990), "A Simple Method for Monotonic Interpolation in One Dimension", http://adsabs.harvard.edu/abs/1990A%26A...239..443S
+
+    Implementation based on http://bl.ocks.org/niclasmattsson/7bceb05fba6c71c78d507adae3d29417
+=#
+
+export
+    LinearMonotonicInterpolation,
+    FiniteDifferenceMonotonicInterpolation,
+    CardinalMonotonicInterpolation,
+    FritschCarlsonMonotonicInterpolation,
+    FritschButlandMonotonicInterpolation,
+    SteffenMonotonicInterpolation
+
+"""
+    MonotonicInterpolationType
+
+Abstract class for all types of monotonic interpolation.
+"""
+abstract type MonotonicInterpolationType <: InterpolationType end
+
+"""
+    LinearMonotonicInterpolation
+
+Simple linear interpolation.
+"""
+struct LinearMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+"""
+    FiniteDifferenceMonotonicInterpolation
+
+Classic cubic interpolation, no tension parameter.
+Finite difference can overshoot for non-monotonic data.
+"""
+struct FiniteDifferenceMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+"""
+    CardinalMonotonicInterpolation(tension)
+
+Cubic cardinal splines, uses `tension` parameter which must be between [0,1]
+Cubin cardinal splines can overshoot for non-monotonic data
+(increasing tension reduces overshoot).
+"""
+struct CardinalMonotonicInterpolation{TTension<:Number} <: MonotonicInterpolationType
+    tension :: TTension # must be in [0, 1]
+end
+
+"""
+    FritschCarlsonMonotonicInterpolation
+
+Monotonic interpolation based on Fritsch & Carlson (1980),
+"Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.
+
+Tangents are first initialized, then adjusted if they are not monotonic
+can overshoot for non-monotonic data
+"""
+struct FritschCarlsonMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+"""
+    FritschButlandMonotonicInterpolation
+
+Monotonic interpolation based on  Fritsch & Butland (1984),
+"A Method for Constructing Local Monotone Piecewise Cubic Interpolants",
+doi:10.1137/0905021.
+
+Faster than FritschCarlsonMonotonicInterpolation (only requires one pass)
+but somewhat higher apparent "tension".
+"""
+struct FritschButlandMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+"""
+    SteffenMonotonicInterpolation
+
+Monotonic interpolation based on Steffen (1990),
+"A Simple Method for Monotonic Interpolation in One Dimension",
+http://adsabs.harvard.edu/abs/1990A%26A...239..443S
+
+Only one pass, results usually between FritschCarlson and FritschButland.
+"""
+struct SteffenMonotonicInterpolation <: MonotonicInterpolationType
+end
+
+"""
+    MonotonicInterpolation
+
+Monotonic interpolation up to third order represented by type, knots and
+coefficients.
+"""
+struct MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType,
+    TKnots<:AbstractVector{<:Number}, TACoeff <: AbstractArray{TEl,1}} <: AbstractInterpolation{T,1,DimSpec{TInterpolationType}}
+
+    it::TInterpolationType
+    knots::TKnots
+    A::TACoeff # constant parts of piecewise polynomials
+    m::Vector{TCoeffs} # coefficients of linear parts of piecewise polynomials
+    c::Vector{TCoeffs} # coefficients of quadratic parts of piecewise polynomials
+    d::Vector{TCoeffs} # coefficients of cubic parts of piecewise polynomials
+end
+
+
+size(A::MonotonicInterpolation) = size(A.knots)
+axes(A::MonotonicInterpolation) = axes(A.knots)
+
+function MonotonicInterpolation(::Type{TWeights}, it::TInterpolationType, knots::TKnots, A::AbstractArray{TEl,1},
+    m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType, TKnots<:AbstractVector{<:Number}}
+
+    isconcretetype(TInterpolationType) || error("The spline type must be a leaf type (was $TInterpolationType)")
+    isconcretetype(TCoeffs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")
+
+    check_monotonic(knots, A)
+
+    cZero = zero(TWeights)
+    if isempty(A)
+        T = Base.promote_op(*, typeof(cZero), eltype(A))
+    else
+        T = typeof(cZero * first(A))
+    end
+
+    MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType, TKnots, typeof(A)}(it, knots, A, m, c, d)
+end
+
+function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::TKnots,
+    A::AbstractArray{TEl,1}, it::TInterpolationType) where {TWeights,TCoeffs,TEl,TKnots<:AbstractVector{<:Number},TInterpolationType<:MonotonicInterpolationType}
+
+    check_monotonic(knots, A)
+
+    # first we need to determine tangents (m)
+    n = length(knots)
+    m, Δ = calcTangents(TCoeffs, knots, A, it)
+    c = Vector{TCoeffs}(undef, n-1)
+    d = Vector{TCoeffs}(undef, n-1)
+    for k ∈ 1:n-1
+        if TInterpolationType == LinearMonotonicInterpolation
+            c[k] = d[k] = zero(TCoeffs)
+        else
+            xdiff = knots[k+1] - knots[k]
+            c[k] = (3*Δ[k] - 2*m[k] - m[k+1]) / xdiff
+            d[k] = (m[k] + m[k+1] - 2*Δ[k]) / (xdiff * xdiff)
+        end
+    end
+
+    MonotonicInterpolation(TWeights, it, knots, A, m, c, d)
+end
+
+function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{TEl,1},
+    it::TInterpolationType) where {TEl,TInterpolationType<:MonotonicInterpolationType}
+
+    interpolate(tweight(A), tcoef(A), knots, A, it)
+end
+
+function (itp::MonotonicInterpolation)(x::Number)
+    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
+    k = searchsortedfirst(itp.knots, x)
+    if k > 1
+        k -= 1
+    end
+    xdiff = x - itp.knots[k]
+    return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
+end
+
+function gradient(itp::MonotonicInterpolation, x::Number)
+    return SVector(gradient1(itp, x))
+end
+
+function gradient1(itp::MonotonicInterpolation, x::Number)
+    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
+    k = searchsortedfirst(itp.knots, x)
+    if k > 1
+        k -= 1
+    end
+    xdiff = x - itp.knots[k]
+    return itp.m[k] + 2*itp.c[k]*xdiff + 3*itp.d[k]*xdiff*xdiff
+end
+
+function hessian(itp::MonotonicInterpolation, x::Number)
+    return SVector(hessian1(itp, x))
+end
+
+function hessian1(itp::MonotonicInterpolation, x::Number)
+    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
+    k = searchsortedfirst(itp.knots, x)
+    if k > 1
+        k -= 1
+    end
+    xdiff = x - itp.knots[k]
+    return 2*itp.c[k] + 6*itp.d[k]*xdiff
+end
+
+@inline function check_monotonic(knots, A)
+    axes(knots) == axes(A) || throw(DimensionMismatch("knot vector must have the same axes as the corresponding array"))
+    issorted(knots) || error("knot-vector must be sorted in increasing order")
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{TEl}, method::LinearMonotonicInterpolation) where {TCoeffs, TEl}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        m[k] = Δk
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{TEl}, method::FiniteDifferenceMonotonicInterpolation) where {TCoeffs, TEl}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        else
+            m[k] = (Δ[k-1] + Δk) / 2
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{TEl}, method::CardinalMonotonicInterpolation{TTension}) where {TTension, TCoeffs, TEl}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        else
+            m[k] = (oneunit(TTension) - method.tension) * (y[k+1] - y[k-1]) / (x[k+1] - x[k-1])
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{TEl}, method::FritschCarlsonMonotonicInterpolation) where {TCoeffs, TEl}
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        else
+            # If any consecutive secant lines change sign (i.e. curve changes direction), initialize the tangent to zero.
+            # This is needed to make the interpolation monotonic. Otherwise set tangent to the average of the secants.
+            if Δk <= zero(Δk)
+                m[k] = zero(TCoeffs)
+            else
+                m[k] = Δ[k-1] * (Δ[k-1] + Δk) / 2.0
+            end
+        end
+    end
+    m[n] = Δ[n-1]
+    #=
+    Fritsch & Carlson derived necessary and sufficient conditions for monotonicity in their 1980 paper. Splines will be
+    monotonic if all tangents are in a certain region of the alpha-beta plane, with alpha and beta as defined below.
+    A robust choice is to put alpha & beta within a circle around origo with radius 3. The FritschCarlson algorithm
+    makes simple initial estimates of tangents and then does another pass over data points to move any outlier tangents
+    into the monotonic region. FritschButland & Steffen algorithms make more elaborate first estimates of tangents that
+    are guaranteed to lie in the monotonic region, so no second pass is necessary.
+    =#
+
+    # Second pass of FritschCarlson: adjust any non-monotonic tangents.
+    for k in 1:n-1
+        Δk = Δ[k]
+        if Δk == zero(TCoeffs)
+            m[k] = zero(TCoeffs)
+            m[k+1] = zero(TCoeffs)
+            continue
+        end
+        α = m[k] / Δk
+        β = m[k+1] / Δk
+        τ = 3.0 / sqrt(α^2 + β^2)
+        if τ < 1.0 # if we're outside the circle with radius 3 then move onto the circle
+            m[k] = τ * α * Δk
+            m[k+1] = τ * β * Δk
+        end
+    end
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y :: AbstractVector{TEl}, method :: FritschButlandMonotonicInterpolation) where {TCoeffs, TEl}
+
+    # based on Fritsch & Butland (1984),
+    # "A Method for Constructing Local Monotone Piecewise Cubic Interpolants",
+    # doi:10.1137/0905021.
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        elseif Δ[k-1] * Δk <= zero(TCoeffs)
+            m[k] = zero(TCoeffs)
+        else
+            α = (1.0 + (x[k+1] - x[k]) / (x[k+1] - x[k-1])) / 3.0
+            m[k] = Δ[k-1] * Δk / (α*Δk + (1.0 - α)*Δ[k-1])
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
+    y::AbstractVector{TEl}, method::SteffenMonotonicInterpolation) where {TCoeffs, TEl}
+
+    # Steffen (1990),
+    # "A Simple Method for Monotonic Interpolation in One Dimension",
+    # http://adsabs.harvard.edu/abs/1990A%26A...239..443S
+
+    n = length(x)
+    Δ = Vector{TCoeffs}(undef, n-1)
+    m = Vector{TCoeffs}(undef, n)
+    for k in 1:n-1
+        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
+        Δ[k] = Δk
+        if k == 1   # left endpoint
+            m[k] = Δk
+        else
+            p = ((x[k+1] - x[k]) * Δ[k-1] + (x[k] - x[k-1]) * Δk) / (x[k+1] - x[k-1])
+            m[k] = (sign(Δ[k-1]) + sign(Δk)) *
+                min(abs(Δ[k-1]), abs(Δk), 0.5*abs(p))
+        end
+    end
+    m[n] = Δ[n-1]
+    return (m, Δ)
+end
+
+# How many non-NoInterp dimensions are there?
+count_interp_dims(::Type{<:MonotonicInterpolationType}) = 1
+
+lbounds(itp::MonotonicInterpolation) = (first(itp.knots),)
+ubounds(itp::MonotonicInterpolation) = (last(itp.knots),)
diff --git a/test/gradient.jl b/test/gradient.jl
index 3e8bd4c5..6c0ec136 100644
--- a/test/gradient.jl
+++ b/test/gradient.jl
@@ -172,4 +172,31 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
     dest = Vector{Float64}(undef, 1)
     Interpolations.gradient!(dest, interp_linear, 1.5)
     @test dest == g
+
+    @testset "Monotonic" begin
+        x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
+        ys = [[-3.0, 0.0, 5.0, 10.0, 18.0, 22.0],
+              [10.0, 0.0, -5.0, 10.0, -8.0, -2.0]]
+        grid = 0.0:0.1:1.0
+
+        itypes = [LinearMonotonicInterpolation(),
+            FiniteDifferenceMonotonicInterpolation(),
+            CardinalMonotonicInterpolation(0.0),
+            CardinalMonotonicInterpolation(0.5),
+            CardinalMonotonicInterpolation(1.0),
+            FritschCarlsonMonotonicInterpolation(),
+            FritschButlandMonotonicInterpolation(),
+            SteffenMonotonicInterpolation()]
+
+        for y in ys
+            for it in itypes
+                itp = interpolate(x, y, it)
+                for t in grid
+                    gradval = epsilon(itp(dual(t, 1.0)))
+                    @test Interpolations.gradient1(itp, t) ≈ gradval atol = 1.e-12
+                    @test Interpolations.gradient(itp, t)[1] ≈ gradval atol = 1.e-12
+                end
+            end
+        end
+    end
 end
diff --git a/test/hessian.jl b/test/hessian.jl
index 2af7574d..297f6965 100644
--- a/test/hessian.jl
+++ b/test/hessian.jl
@@ -1,4 +1,4 @@
-using Test, Interpolations, LinearAlgebra
+using Test, Interpolations, LinearAlgebra, ForwardDiff
 
 @testset "Hessians" begin
     nx = 5
@@ -42,4 +42,32 @@ using Test, Interpolations, LinearAlgebra
     v = A2[:, 2]
     itpcol = interpolate(v, BSpline(Quadratic(Flat(OnCell()))))
     @test Interpolations.hessian(itp, 3.2, 2) == Interpolations.hessian(itpcol, 3.2)
+
+
+    @testset "Monotonic" begin
+        x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
+        ys = [[-3.0, 0.0, 5.0, 10.0, 18.0, 22.0],
+              [10.0, 0.0, -5.0, 10.0, -8.0, -2.0]]
+        grid = 0.0:0.1:1.0
+
+        itypes = [LinearMonotonicInterpolation(),
+            FiniteDifferenceMonotonicInterpolation(),
+            CardinalMonotonicInterpolation(0.0),
+            CardinalMonotonicInterpolation(0.5),
+            CardinalMonotonicInterpolation(1.0),
+            FritschCarlsonMonotonicInterpolation(),
+            FritschButlandMonotonicInterpolation(),
+            SteffenMonotonicInterpolation()]
+
+        for y in ys
+            for it in itypes
+                itp = interpolate(x, y, it)
+                for t in grid
+                    hessval = ForwardDiff.hessian(u -> itp(u[1]), [t])[1, 1]
+                    @test Interpolations.hessian1(itp, t) ≈ hessval atol = 1.e-12
+                    @test Interpolations.hessian(itp, t)[1] ≈ hessval atol = 1.e-12
+                end
+            end
+        end
+    end
 end
diff --git a/test/io.jl b/test/io.jl
index 6027c13e..5c5a9768 100644
--- a/test/io.jl
+++ b/test/io.jl
@@ -57,6 +57,13 @@ using Test
         @test summary(sitp2) == "21×41 scale(interpolate(OffsetArray(::Array{Float64,2}, 0:22, 0:42), BSpline(Quadratic(Flat(OnGrid())))), (-5.0:0.5:5.0,$SPACE-4.0:0.2:4.0)) with element type Float64"
     end
 
+    @testset "Monotonic" begin
+        A = rand(20)
+        A_x = collect(1.0:2.0:40.0)
+        itp = interpolate(A_x, A, FritschButlandMonotonicInterpolation())
+        @test summary(itp) == "20-element interpolate(::Array{Float64,1}, ::Array{Float64,1}, FritschButlandMonotonicInterpolation()) with element type Float64"
+    end
+
     @testset "Extrapolation" begin
         A = rand(8,20)
 
@@ -67,4 +74,4 @@ using Test
         etpf = extrapolate(itpg, NaN)
         @test summary(etpf) == "8×20 extrapolate(interpolate(::Array{Float64,2}, BSpline(Linear())), NaN) with element type Float64"
     end
-end # Module
+end # testset
diff --git a/test/monotonic/runtests.jl b/test/monotonic/runtests.jl
new file mode 100644
index 00000000..9d71d221
--- /dev/null
+++ b/test/monotonic/runtests.jl
@@ -0,0 +1,54 @@
+
+using Test
+using Interpolations
+
+@testset "Monotonic" begin
+    x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
+    y = [(true, [-3.0, 0.0, 5.0, 10.0, 18.0, 22.0]),
+         (false, [10.0, 0.0, -5.0, 10.0, -8.0, -2.0])]
+
+    # second item indicates if interpolating function can overshoot
+    # for non-monotonic data
+    itypes = [(LinearMonotonicInterpolation(), false),
+        (FiniteDifferenceMonotonicInterpolation(), true),
+        (CardinalMonotonicInterpolation(0.0), true),
+        (CardinalMonotonicInterpolation(0.5), true),
+        (CardinalMonotonicInterpolation(1.0), false),
+        (FritschCarlsonMonotonicInterpolation(), true),
+        (FritschButlandMonotonicInterpolation(), false),
+        (SteffenMonotonicInterpolation(), false)]
+
+    for (it, overshoot) in itypes
+        for yi in 1:length(y)
+            monotonic, ys = y[yi]
+            itp = interpolate(x, ys, it)
+            for j in 1:6
+                # checking values at nodes
+                @test itp(x[j]) ≈ ys[j] rtol = 1.e-15 atol = 1.e-14
+
+                # checking overshoot for non-monotonic data
+                # and monotonicity for monotonic data
+                if !monotonic && overshoot
+                    continue
+                end
+                if j < 6
+                    r = range(x[j], stop = x[j+1], length = 100)
+                    for k in 1:length(r)-1
+                        @test (ys[j+1] - ys[j]) * (itp(r[k+1]) - itp(r[k])) >= 0.0
+                    end
+                end
+            end
+        end
+    end
+
+    # fail tests
+    xWrong = [0.0, 1.0, -5.0, 3.0]
+    y2 = [1.0, 2.0, 3.0, 4.0]
+    y3 = [-3.0 0.0 5.0 10.0 18.0 22.0]
+
+    for (it, overshoot) in itypes
+        @test_throws ErrorException interpolate(xWrong, y2, it)
+        @test_throws DimensionMismatch interpolate(x, y2, it)
+        @test_throws MethodError interpolate(x, y3, it)
+    end
+end
diff --git a/test/plots/monotonic/sanity.jl b/test/plots/monotonic/sanity.jl
new file mode 100644
index 00000000..419e84c1
--- /dev/null
+++ b/test/plots/monotonic/sanity.jl
@@ -0,0 +1,25 @@
+using Plots
+using Interpolations
+
+pyplot()
+
+x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
+y = [10.0, 0.0, -5.0, 10.0, -8.0, -2.0]
+grid = 0.0:0.002:1.0
+
+itypes = [LinearMonotonicInterpolation(),
+    FiniteDifferenceMonotonicInterpolation(),
+    CardinalMonotonicInterpolation(0.0),
+    CardinalMonotonicInterpolation(0.5),
+    CardinalMonotonicInterpolation(1.0),
+    FritschCarlsonMonotonicInterpolation(),
+    FritschButlandMonotonicInterpolation(),
+    SteffenMonotonicInterpolation()
+    ]
+
+scatter(x, y, label = "data")
+for itype in itypes
+    itp = interpolate(x, y, itype)
+    plot!([x for x in grid], [itp(x) for x in grid], label = string(itype))
+end
+gui()
diff --git a/test/runtests.jl b/test/runtests.jl
index c31533d9..8d37d5ac 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -25,6 +25,9 @@ const isci = get(ENV, "CI", "") in ("true", "True")
     # scaling tests
     include("scaling/runtests.jl")
 
+    # monotonic tests
+    include("monotonic/runtests.jl")
+
     # test gradient evaluation
     include("gradient.jl")
     isci && println("finished gradient")