diff --git a/Project.toml b/Project.toml
index e135ba06..dcc05fc9 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.18.0"
+version = "0.18.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/docs/src/api.md b/docs/src/api.md
index c33907b4..cd40cda2 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -47,6 +47,7 @@ constant_term
 linear_polynomial
 nonlinear_polynomial
 inverse
+inverse_map
 abs
 norm
 isapprox
diff --git a/docs/src/index.md b/docs/src/index.md
index d7ac947c..6f4bef5b 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -49,6 +49,7 @@ programs in Physics and in Mathematics at UNAM, during the second half of 2013.
 We thank the participants of the course for putting up with the half-baked
 material and contributing energy and ideas.
 
-We acknowledge financial support from DGAPA-UNAM PAPIME grants PE-105911 and
-PE-107114, and PAPIIT grants IG-101113 and IG-100616. LB acknowledges
-support through a *Cátedra Marcos Moshinsky* (2013).
+We acknowledge financial support from DGAPA-UNAM PAPIME grants
+PE-105911 and PE-107114, and DGAPA-PAPIIT grants IG-101113,
+IG-100616, IG-100819 and IG-101122.
+LB acknowledges support through a *Cátedra Marcos Moshinsky* (2013).
diff --git a/docs/src/userguide.md b/docs/src/userguide.md
index 30446a07..184a51be 100644
--- a/docs/src/userguide.md
+++ b/docs/src/userguide.md
@@ -193,6 +193,8 @@ differentiate(5, exp(shift_taylor(1.0))) == exp(1.0)    # 5-th differentiate of
 derivative(exp(1+t), 3)    # Taylor1 polynomial of the 3-rd derivative of `exp(1+t)`
 ```
 
+We also have methods to invert a `Taylor1` polynomial, using either [`inverse`](@ref), which uses Lagrange inversion, or [`inverse_map`](@ref), which uses an algorithm developed by M. Berz; the latter can also be used for `TaylorN` polynomials; see below.
+
 To evaluate a Taylor series at a given point, Horner's rule is used via the
 function `evaluate(a, dt)`. Here, `dt` is the increment from
 the point ``t_0`` around which the Taylor expansion of `a` is calculated,
@@ -230,7 +232,7 @@ The former returns
 the expansion of a function around a given value `t0`, mimicking the use
 of `shift_taylor` above. In turn, `update!`
 provides an in-place update of a given Taylor polynomial, that is, it shifts
-it further by the provided amount.
+it further by the provided amount. Note that the type of the `Taylor1` polynomial and the shifted point must be compatible, in the sense that the latter must be convertable into the parametric type of the former.
 
 ```@repl userguide
 p = taylor_expand( x -> sin(x), pi/2, order=16) # 16-th order expansion of sin(t) around pi/2
diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 2e508e66..8df08885 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -54,7 +54,7 @@ import Base.float
 export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries, TS
 
 export getcoeff, derivative, integrate, differentiate,
-    evaluate, evaluate!, inverse, set_taylor1_varname,
+    evaluate, evaluate!, inverse, inverse_map, set_taylor1_varname,
     show_params_TaylorN, show_monomials, displayBigO, use_show_default,
     get_order, get_numvars,
     set_variables, get_variables,
diff --git a/src/evaluate.jl b/src/evaluate.jl
index 6594c3fa..a496c607 100644
--- a/src/evaluate.jl
+++ b/src/evaluate.jl
@@ -331,8 +331,7 @@ end
 
 function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number}
     R = promote_type(T, typeof(vals[1]))
-    a_length = length(a)
-    suma = zeros(R, a_length)
+    suma = zeros(R, length(a))
     @inbounds for homPol in eachindex(a)
         suma[homPol+1] = _evaluate(a[homPol], vals)
     end
@@ -499,18 +498,20 @@ end
 
 function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1},
         x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN}
-    x0 .= _evaluate.( x, Ref(δx), Ref(Val(sorting)) )
+    x0 .= evaluate.( x, Ref(δx), sorting = sorting)
     return nothing
 end
 
-function evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T}, valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
+function evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T},
+        valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
     @inbounds for homPol in eachindex(a)
         _evaluate!(dest, a[homPol], vals, valscache, aux)
     end
     return nothing
 end
 
-function evaluate!(a::AbstractArray{TaylorN{T}}, vals::NTuple{N,TaylorN{T}}, dest::AbstractArray{TaylorN{T}}) where {N,T<:Number}
+function evaluate!(a::AbstractArray{TaylorN{T}}, vals::NTuple{N,TaylorN{T}},
+        dest::AbstractArray{TaylorN{T}}) where {N,T<:Number}
     # initialize evaluation cache
     valscache = [zero(val) for val in vals]
     aux = zero(dest[1])
diff --git a/src/functions.jl b/src/functions.jl
index bb27b9af..ae54c6d6 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -1228,7 +1228,7 @@ end
     inverse(f)
 
 Return the Taylor expansion of ``f^{-1}(t)``, of order `N = f.order`,
-for `f::Taylor1` polynomial if the first coefficient of `f` is zero.
+for `f::Taylor1` polynomial, assuming the first coefficient of `f` is zero.
 Otherwise, a `DomainError` is thrown.
 
 The algorithm implements Lagrange inversion at ``t=0`` if ``f(0)=0``:
@@ -1246,24 +1246,85 @@ function inverse(f::Taylor1{T}) where {T<:Number}
     if !iszero(f[0])
         throw(DomainError(f,
         """
-        Evaluation of Taylor1 series at 0 is non-zero. For high accuracy, revert
-        a Taylor1 series with first coefficient 0 and re-expand about f(0).
+        Evaluation of Taylor1 series at 0 is non-zero; revert
+        a Taylor1 series with constant coefficient 0 and re-expand about f(0).
         """))
     end
-    z = Taylor1(T,f.order)
+    z = Taylor1(T, f.order)
     zdivf = z/f
     zdivfpown = zdivf
-    S = TS.numtype(zdivf)
-    coeffs = zeros(S,f.order+1)
+    res = Taylor1(zero(TS.numtype(zdivf)), f.order)
 
-    @inbounds for n in 1:f.order
-        coeffs[n+1] = zdivfpown[n-1]/n
+    @inbounds for ord in 1:f.order
+        res[ord] = zdivfpown[ord-1]/ord
         zdivfpown *= zdivf
     end
-    Taylor1(coeffs, f.order)
+    return res
 end
 
 
+@doc doc"""
+    inverse_map(f)
+
+Return the Taylor expansion of ``f^{-1}(t)``, of order `N = f.order`,
+for `Taylor1` or `TaylorN` polynomials, assuming the first coefficient of `f` is zero.
+Otherwise, a `DomainError` is thrown.
+
+This method is based in the algorithm by M. Berz, Modern map methods in
+Particle Beam Physics, Academic Press (1999), Sect 2.3.1.
+See [`inverse`](@ref) (for `f::Taylor1`).
+
+"""
+function inverse_map(p::Taylor1)
+    if !iszero(constant_term(p))
+        throw(DomainError(p,
+        """
+        Evaluation of Taylor1 series at 0 is non-zero; revert
+        a Taylor1 series with constant coefficient 0 and re-expand about f(0).
+        """))
+    end
+    inv_m_pol = inv(linear_polynomial(p)[1])
+    n_pol = inv_m_pol * nonlinear_polynomial(p)
+    scaled_ident = inv_m_pol * Taylor1(p.order)
+    res = scaled_ident
+    aux1 = zero(res)
+    aux2 = zero(res)
+    for ord in 1:p.order
+        _horner!(aux2, n_pol, res, aux1)
+        subst!(res, scaled_ident, aux2, ord)
+    end
+    return res
+end
+
+function inverse_map(p::Vector{TaylorN{T}}) where {T<:NumberNotSeries}
+    if !iszero(constant_term(p))
+        throw(DomainError(p,
+        """
+        Evaluation of Taylor1 series at 0 is non-zero; revert
+        a Taylor1 series with constant coefficient 0 and re-expand about f(0).
+        """))
+    end
+    @assert length(p) == get_numvars()
+    inv_m_pol = inv(jacobian(p))
+    n_pol = inv_m_pol * nonlinear_polynomial(p)
+    scaled_ident = inv_m_pol * TaylorN.(1:get_numvars(), order=get_order(p[1]))
+    res = deepcopy(scaled_ident)
+    aux = zero.(res)
+    auxvec = [zero(res[1]) for val in eachindex(res[1])]
+    valscache = [zero(val) for val in res]
+    aaux = zero(res[1])
+    for ord in 1:get_order(p[1])
+        t_res = (res...,)
+        for i = 1:get_numvars()
+            zero!.(auxvec)
+            _evaluate!(auxvec, n_pol[i], t_res, valscache, aaux)
+            aux[i] = sum( auxvec )
+            subst!(res[i], scaled_ident[i], aux[i], ord)
+        end
+    end
+    return res
+end
+
 
 # Documentation for the recursion relations
 @doc doc"""
diff --git a/src/other_functions.jl b/src/other_functions.jl
index 8d225851..3d7e6be9 100644
--- a/src/other_functions.jl
+++ b/src/other_functions.jl
@@ -254,15 +254,23 @@ end
 
 Takes `a <: Union{Taylo1,TaylorN}` and expands it around the coordinate `x0`.
 """
-function update!(a::Taylor1, x0::T) where {T<:Number}
+function update!(a::Taylor1{T}, x0::T) where {T<:Number}
     a.coeffs .= evaluate(a, Taylor1([x0, one(x0)], a.order) ).coeffs
-    nothing
+    return nothing
+end
+function update!(a::Taylor1{T}, x0::S) where {T<:Number, S<:Number}
+    xx0 = convert(T, x0)
+    return update!(a, xx0)
 end
 
 #update! function for TaylorN
-function update!(a::TaylorN, vals::Vector{T}) where {T<:Number}
+function update!(a::TaylorN{T}, vals::Vector{T}) where {T<:Number}
     a.coeffs .= evaluate(a, get_variables(a.order) .+ vals).coeffs
-    nothing
+    return nothing
+end
+function update!(a::TaylorN{T}, vals::Vector{S}) where {T<:Number, S<:Number}
+    vv = convert(Vector{T}, vals)
+    return update!(a, vv)
 end
 
 function update!(a::Union{Taylor1,TaylorN})
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index be842a61..8e49a1c0 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -344,6 +344,20 @@ end
     txy[2:end-1] .= ( 1.0 - xT*yT + 0.5*xT^2*yT - (2/3)*xT*yT^3 - 0.5*xT^2*yT^2  + 7*xT^3*yT )[2:end-1]
     @test txy[2:end-1] == ( 1.0 - xT*yT + 0.5*xT^2*yT - (2/3)*xT*yT^3 - 0.5*xT^2*yT^2  + 7*xT^3*yT )[2:end-1]
 
+    ident = [xT, yT]
+    pN = [x+y, x-y]
+    @test evaluate.(inverse_map(pN), Ref(pN)) == ident
+    @test evaluate.(pN, Ref(inverse_map(pN))) == ident
+    pN = [exp(xT)-1, log(1+yT)]
+    @test inverse_map(pN) ≈ [log(1+xT), exp(yT)-1]
+    @test evaluate.(pN, Ref(inverse_map(pN))) ≈ ident
+    pN = [tan(xT), atan(yT)]
+    @test evaluate.(inverse_map(pN), Ref(pN)) ≈ ident
+    @test evaluate.(pN, Ref(inverse_map(pN))) ≈ ident
+    pN = [sin(xT), asin(yT)]
+    @test evaluate.(inverse_map(pN), Ref(pN)) ≈ ident
+    @test evaluate.(pN, Ref(inverse_map(pN))) ≈ ident
+
     a = -5.0 + sin(xT+yT^2)
     b = deepcopy(a)
     @test a[:] == a[0:end]
@@ -749,7 +763,7 @@ end
     xysq = x^2 + y^2
     update!(xysq,[1.0,-2.0])
     @test xysq == (x+1.0)^2 + (y-2.0)^2
-    update!(xysq,[-1.0,2.0])
+    update!(xysq,[-1,2])
     @test xysq == x^2 + y^2
 
     #test function-like behavior for TaylorN
diff --git a/test/onevariable.jl b/test/onevariable.jl
index 847ff8b3..08b195b1 100644
--- a/test/onevariable.jl
+++ b/test/onevariable.jl
@@ -259,7 +259,7 @@ Base.iszero(::SymbNumber) = false
     tsq = t^2
     update!(tsq,2.0)
     @test tsq == (t+2.0)^2
-    update!(tsq,-2.0)
+    update!(tsq,-2)
     @test tsq == t^2
 
     @test log(exp(tsquare)) == tsquare
@@ -514,9 +514,14 @@ Base.iszero(::SymbNumber) = false
 
     @test promote(ta(0.0), t) == (ta(0.0),ta(0.0))
 
-    @test norm((inverse(exp(t)-1) - log(1+t)).coeffs) < 2tol1
+    @test inverse(exp(t)-1) ≈ log(1+t)
     cfs = [(-n)^(n-1)/factorial(n) for n = 1:15]
     @test norm(inverse(t*exp(t))[1:end]./cfs .- 1) < 4tol1
+    @test inverse(tan(t))(tan(t)) ≈ t
+    @test atan(inverse(atan(t))) ≈ t
+    @test inverse_map(sin(t))(sin(t)) ≈ t
+    @test sinh(inverse_map(sinh(t))) ≈ t
+    @test inverse_map(tanh(t)) ≈ inverse(tanh(t))
 
     @test_throws ArgumentError Taylor1([1,2,3], -2)
     @test_throws DomainError abs(ta(big(0)))