diff --git a/Project.toml b/Project.toml
index ffbcaffd..9bf29e50 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.19.1"
+version = "0.20.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index 8d06dfb4..bc29b0ba 100644
--- a/ext/TaylorSeriesIAExt.jl
+++ b/ext/TaylorSeriesIAExt.jl
@@ -4,7 +4,8 @@ using TaylorSeries
 
 import Base: ^, sqrt, log, asin, acos, acosh, atanh, iszero, ==
 
-import TaylorSeries: _pow, pow!, evaluate, _evaluate, normalize_taylor, aff_normalize
+import TaylorSeries: _pow, pow!, evaluate, _evaluate, _evaluate!,
+        normalize_taylor, aff_normalize
 
 using IntervalArithmetic
 
@@ -751,7 +752,12 @@ for T in (:Taylor1, :TaylorN)
 end
 
 
-function evaluate(a::Taylor1{T}, dx::Interval{S}) where {T<:Real, S<:NumTypes}
+function evaluate(a::Taylor1{Interval{T}}, dx::S) where {T<:NumTypes, S<:NumTypes}
+    dxI, _ = promote(dx, constant_term(a))
+    return evaluate(a, dxI)
+end
+
+function evaluate(a::Taylor1{T}, dx::Interval{T}) where {T<:NumTypes}
     order = a.order
     uno = one(dx)
     dx2 = dx^2
@@ -771,16 +777,44 @@ function evaluate(a::Taylor1{T}, dx::Interval{S}) where {T<:Real, S<:NumTypes}
     return sum_even + sum_odd*dx
 end
 
-function evaluate(a::TaylorN{T}, dx::AbstractVector{Interval{S}}) where
-        {T<:Real, S<:NumTypes}
-    @assert length(dx) == get_numvars()
-    suma = zero(constant_term(a)) + interval(zero(T))
-    @inbounds for homPol in reverse(eachindex(a))
-        suma += evaluate(a.coeffs[homPol+1], dx)
+function evaluate(a::Taylor1{Interval{T}}, dx::Interval{T}) where {T<:NumTypes}
+    order = a.order
+    uno = one(dx)
+    dx2 = dx^2
+    if iseven(order)
+        kend = order-2
+        @inbounds sum_even = a[end]*uno
+        @inbounds sum_odd = a[end-1]*zero(dx)
+    else
+        kend = order-3
+        @inbounds sum_odd = a[end]*uno
+        @inbounds sum_even = a[end-1]*uno
     end
-    return suma
+    @inbounds for k in kend:-2:0
+        sum_odd = sum_odd*dx2 + a[k+1]
+        sum_even = sum_even*dx2 + a[k]
+    end
+    return sum_even + sum_odd*dx
+end
+
+function evaluate(a::TaylorN{Interval{T}}, dx::AbstractVector{S}) where {T<:NumTypes, S<:NumTypes}
+    @assert length(dx) == get_numvars()
+    return evaluate(a, (dx...,); sorting=false)
 end
 
+function evaluate(a::TaylorN{T}, dx::AbstractVector{Interval{T}}) where {T<:NumTypes}
+    @assert length(dx) == get_numvars()
+    return evaluate(a, (dx...,); sorting=false)
+end
+
+function evaluate(a::TaylorN{Interval{T}}, dx::AbstractVector{Interval{T}}) where {T<:NumTypes}
+    @assert length(dx) == get_numvars()
+    return evaluate(a, (dx...,); sorting=false)
+end
+
+evaluate(a::TaylorN{T}, vals::AbstractVector{TaylorN{Interval{T}}}) where
+        {T<:NumTypes} = evaluate(a, (vals...,); sorting=false)
+
 function evaluate(a::Taylor1{TaylorN{T}}, dx::Interval{S}) where {T<:Real, S<:Real}
     order = a.order
     uno = one(dx)
@@ -874,6 +908,17 @@ function TS._evaluate(a::TaylorN{T}, vals::NTuple{N,TaylorN{Interval{S}}}) where
     return suma
 end
 
+function TS._evaluate(a::TaylorN{Interval{T}}, vals::NTuple{N,TaylorN{Interval{T}}}) where
+        {N, T<:NumTypes}
+    @assert get_numvars() == N
+    a_length = length(a)
+    suma = zeros(T, a_length)
+    @inbounds for homPol in 1:a_length
+        suma[homPol] = TS._evaluate(a.coeffs[homPol], vals)
+    end
+    return suma
+end
+
 function TS._evaluate(a::HomogeneousPolynomial{T},
         vals::NTuple{N,TaylorN{Interval{S}}}) where {N, T<:Real, S<:NumTypes}
     ct = TS.coeff_table[a.order+1]
@@ -907,80 +952,30 @@ are mapped by an affine transformation to the intervals `-1..1`
 normalize_taylor(a::TaylorN, I::AbstractVector{Interval{T}},
     symI::Bool=true) where {T<:NumTypes} = _normalize(a, I, Val(symI))
 
-aff_normalize(x, I::Interval, ::Val{true})  = interval(mid(I)) + x*interval(radius(I))
-aff_normalize(x, I::Interval, ::Val{false}) = interval(inf(I)) + x*interval(diam(I))
-
-#  I -> -1..1
-function _normalize(a::Taylor1, I::Interval{T}, ::Val{true}) where {T<:NumTypes}
-    t = Taylor1(TS.numtype(a), get_order(a))
-    return a(aff_normalize(t, I, Val(true)))
-end
+aff_normalize(x, I::Interval, ::Val{true})  = interval(mid(I)) + x * interval(radius(I))
+aff_normalize(x, I::Interval, ::Val{false}) = interval(inf(I)) + x * interval(diam(I))
 
-# function _normalize(a::Taylor1{Interval{T}}, I::Interval{T}, ::Val{true}) where
-#         {T<:NumTypes}
-#     t = Taylor1(TS.numtype(a), get_order(a))
-#     return a(aff_normalize(t, I, Val(true)))
-# end
-
-#  I -> 0..1
-function _normalize(a::Taylor1, I::Interval{T}, ::Val{false}) where {T<:NumTypes}
-    t = Taylor1(TS.numtype(a), get_order(a))
-    return a(aff_normalize(t, I, Val(false)))
-end
-
-# function _normalize(a::Taylor1{Interval{T}}, I::Interval{T}, ::Val{false}) where
-#         {T<:NumTypes}
-#     t = Taylor1(TS.numtype(a), get_order(a))
-#     return a(aff_normalize(t, I, Val(false)))
-# end
-
-#  I -> IntervalBox(-1..1, Val(N))
-function _normalize(a::TaylorN, I::AbstractVector{Interval{T}},
-        ::Val{true}) where {T<:NumTypes}
-    order = get_order(a)
-    x = Vector{typeof(a)}(undef, length(I))
-    @inbounds for ind in eachindex(x)
-        # x[ind] = mid(I[ind]) + TaylorN(ind, order=order)*radius(I[ind])
-        x[ind] = aff_normalize(TaylorN(ind, order=order), I[ind], Val(true))
+for bb in (:true, :false)
+    @eval function _normalize(a::Taylor1, I::Interval{T}, ::Val{$bb}) where {T<:NumTypes}
+        S = promote_type(TS.numtype(a), Interval{T})
+        t = Taylor1(S, get_order(a))
+        return a(aff_normalize(t, I, Val($bb)))
     end
-    return a(x)
-end
 
-# function _normalize(a::TaylorN{Interval{T}}, I::AbstractVector{Interval{T}},
-#         ::Val{true}) where {T<:NumTypes}
-#     order = get_order(a)
-#     x = Vector{typeof(a)}(undef, length(I))
-#     @inbounds for ind in eachindex(x)
-#         # x[ind] = interval(mid(I[ind])) +
-#         #     TaylorN(ind, order=order)*interval(radius(I[ind]))
-#         x[ind] = aff_normalize(TaylorN(ind, order=order), I[ind], Val(true))
-#     end
-#     return a(x)
-# end
-
-#  I -> IntervalBox(0..1, Val(N))
-function _normalize(a::TaylorN, I::AbstractVector{Interval{T}},
-        ::Val{false}) where {T<:NumTypes}
-    order = get_order(a)
-    x = Vector{typeof(a)}(undef, length(I))
-    @inbounds for ind in eachindex(x)
-        # x[ind] = inf(I[ind]) + TaylorN(ind, order=order)*diam(I[ind])
-        x[ind] = aff_normalize(TaylorN(ind, order=order), I[ind], Val(false))
+    @eval function _normalize(a::TaylorN, I::AbstractVector{Interval{T}},
+            ::Val{$bb}) where {T<:NumTypes}
+        order = get_order(a)
+        S = promote_type(TS.numtype(a), Interval{T})
+        x = Vector{TaylorN{S}}(undef, length(I))
+        @inbounds for ind in eachindex(x)
+            # x[ind] = mid(I[ind]) + TaylorN(ind, order=order)*radius(I[ind])
+            x[ind] = aff_normalize(TaylorN(S, ind, order=order), I[ind], Val($bb))
+        end
+        aa = convert(TaylorN{S}, a)
+        return evaluate(aa, x)
     end
-    return a(x)
 end
 
-# function _normalize(a::TaylorN{Interval{T}}, I::AbstractVector{Interval{T}},
-#         ::Val{false}) where {T<:NumTypes}
-#     order = get_order(a)
-#     x = Vector{typeof(a)}(undef, length(I))
-#     @inbounds for ind in eachindex(x)
-#         # x[ind] = interval(inf(I[ind])) +
-#         #     TaylorN(ind, order=order)*interval(diam(I[ind]))
-#         x[ind] = aff_normalize(TaylorN(ind, order=order), I[ind], Val(false))
-#     end
-#     return a(x)
-# end
 
 # Printing-related methods numbr2str
 function TS.numbr2str(zz::Interval, ifirst::Bool=false)
diff --git a/src/evaluate.jl b/src/evaluate.jl
index 4a8df6fb..3ef31fd9 100644
--- a/src/evaluate.jl
+++ b/src/evaluate.jl
@@ -277,13 +277,13 @@ function evaluate(a::TaylorN, vals::NTuple{N,<:AbstractSeries};
 end
 
 evaluate(a::TaylorN{T}, vals::AbstractVector{<:Number}; sorting::Bool=true) where
-    {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting=sorting)
+    {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting)
 
 evaluate(a::TaylorN{T}, vals::AbstractVector{<:AbstractSeries}; sorting::Bool=false) where
-    {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting=sorting)
+    {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting)
 
 evaluate(a::TaylorN{Taylor1{T}}, vals::AbstractVector{S};
-    sorting::Bool=false) where {T, S} = evaluate(a, (vals...,); sorting=sorting)
+    sorting::Bool=false) where {T, S} = evaluate(a, (vals...,); sorting)
 
 function evaluate(a::TaylorN{T}, s::Symbol, val::S) where
         {T<:Number, S<:NumberNotSeriesN}
@@ -292,8 +292,8 @@ function evaluate(a::TaylorN{T}, s::Symbol, val::S) where
     return evaluate(a, ind, val)
 end
 
-function evaluate(a::TaylorN{T}, ind::Int, val::S) where
-        {T<:Number, S<:NumberNotSeriesN}
+function evaluate(a::TaylorN{T}, ind::Int, val::S) where {T<:Number,
+        S<:NumberNotSeriesN}
     @assert (1 ≤ ind ≤ get_numvars()) "Invalid `ind`; it must be between 1 and `get_numvars()`"
     R = promote_type(T,S)
     return _evaluate(convert(TaylorN{R}, a), ind, convert(R, val))
@@ -325,8 +325,8 @@ _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where
 _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} =
     sum( _evaluate(a, vals) )
 
-function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN},
-        ::Val{false}) where {N,T<:Number}
+function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}, ::Val{false}) where
+        {N,T<:Number}
     R = promote_type(T, TS.numtype(vals[1]))
     res = TaylorN(zero(R), vals[1].order)
     vvals = ntuple(i -> convert(TaylorN{R}, vals[i]), length(vals))
diff --git a/test/runtests.jl b/test/runtests.jl
index ef447b14..1057249a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,3 +21,6 @@ testfiles = (
 for file in testfiles
     include(file)
 end
+
+# Back to TaylorN default parameters
+set_variables("x", order=6, numvars=2);