diff --git a/docs/src/userguide.md b/docs/src/userguide.md
index d89181a6..8bd85cfc 100644
--- a/docs/src/userguide.md
+++ b/docs/src/userguide.md
@@ -142,6 +142,23 @@ eBig = evaluate( exp(tBig), one(BigFloat) )
 e - eBig
 ```
 
+Another way to obtain the value of a `Taylor1` polynomial `p` at a given value `x`, is to call `p` as if it were a function:
+
+```@repl userguide
+t = Taylor1(15)
+p = sin(t)
+evaluate(p, pi/2) #get value of p at pi/2 using `evaluate`
+p(pi/2) #get value of p at pi/2 by calling p as a function
+p(pi/2) == evaluate(p, pi/2)
+p(0.0) #get value of `p` at 0.0 by calling p as a function
+evaluate(p) #get p 0-th order coefficient using `evaluate`
+p() #a shortcut to get 0-th order coefficient of `p`
+p() == evaluate(p)
+p() == p(0.0)
+```
+
+Note that the syntax `p(x)` is equivalent to `evaluate(p,x)`, whereas `p()` is equivalent to `evaluate(p)`. For more details about function-like behavior for a given type in Julia, see the [Function-like objects](https://docs.julialang.org/en/stable/manual/methods/#Function-like-objects-1) section of the Julia manual.
+
 ## Many variables
 
 A polynomial in ``N>1`` variables can be represented in (at least) two ways:
@@ -283,6 +300,15 @@ number of independent variables.
 evaluate(exy, [.1,.02]) == e^0.12
 ```
 
+Analogously to `Taylor1`, another way to obtain the value of a `TaylorN` polynomial `p` at a given point `x`, is to call `p` as if it were a function:
+
+```@repl userguide
+exy([.1,.02])
+exy([.1,.02]) == e^0.12
+```
+
+Again, note that the syntax `p(x)` for `p::TaylorN` is equivalent to `evaluate(p,x)`, whereas `p()` is equivalent to `evaluate(p)`. For more details about function-like behavior for a given type in Julia, see the [Function-like objects](https://docs.julialang.org/en/stable/manual/methods/#Function-like-objects-1) section of the Julia manual.
+
 Functions to compute the gradient, Jacobian and
 Hessian have also been implemented. Using the
 polynomials ``p = x^3 + 2x^2 y - 7 x + 2`` and ``q = y-x^4`` defined above,
@@ -305,7 +331,7 @@ Other specific applications are described in the next [section](examples).
 ## Mixtures
 
 As mentioned above, `Taylor1{T}`, `HomogeneousPolynomial{T}` and `TaylorN{T}`
-are parameterized structures structures such that `T<:AbstractSeries`, the latter
+are parameterized structures such that `T<:AbstractSeries`, the latter
 is a subtype of `Number`. Then, we may actually define Taylor expansions in
 ``N+1`` variables, where one of the variables (the `Taylor1` variable) is
 somewhat special.
diff --git a/src/constructors.jl b/src/constructors.jl
index f0724815..bd0782ce 100644
--- a/src/constructors.jl
+++ b/src/constructors.jl
@@ -28,6 +28,9 @@ DataType for polynomial expansions in one independent variable.
 - `coeffs :: Array{T,1}` Expansion coefficients; the $i$-th
     component is the coefficient of degree $i-1$ of the expansion.
 - `order  :: Int64` Maximum order (degree) of the polynomial.
+
+Note that `Taylor1` variables are callable. For more information, see
+[`evaluate`](@ref).
 """
 immutable Taylor1{T<:Number} <: AbstractSeries{T}
     coeffs :: Array{T,1}
@@ -81,6 +84,9 @@ DataType for homogenous polynomials in many (>1) independent variables.
 polynomial; the $i$-th component is related to a monomial, where the degrees
 of the independent variables are specified by `coeff_table[order+1][i]`.
 - `order   :: Int` order (degree) of the homogenous polynomial.
+
+Note that `HomogeneousPolynomial` variables are callable. For more information,
+see [`evaluate`](@ref).
 """
 immutable HomogeneousPolynomial{T<:Number} <: AbstractSeries{T}
     coeffs  :: Array{T,1}
@@ -137,6 +143,9 @@ DataType for polynomial expansions in many (>1) independent variables.
 `HomogeneousPolynomial` entries. The $i$-th component corresponds to the
 homogeneous polynomial of degree $i-1$.
 - `order   :: Int`  maximum order of the polynomial expansion.
+
+Note that `TaylorN` variables are callable. For more information, see
+[`evaluate`](@ref).
 """
 immutable TaylorN{T<:Number} <: AbstractSeries{T}
     coeffs  :: Array{HomogeneousPolynomial{T},1}
diff --git a/src/evaluate.jl b/src/evaluate.jl
index adf26561..0203ea14 100644
--- a/src/evaluate.jl
+++ b/src/evaluate.jl
@@ -11,7 +11,9 @@
     evaluate(a, [dx])
 
 Evaluate a `Taylor1` polynomial using Horner's rule (hand coded). If `dx` is
-ommitted, its value is considered as zero.
+ommitted, its value is considered as zero. Note that a `Taylor1` polynomial `a`
+may also be evaluated by calling it as a function; that is, the syntax `a(dx)`
+is equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`.
 """
 function evaluate{T<:Number}(a::Taylor1{T}, dx::T)
     @inbounds suma = a[end]
@@ -20,6 +22,7 @@ function evaluate{T<:Number}(a::Taylor1{T}, dx::T)
     end
     suma
 end
+
 function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, dx::S)
     R = promote_type(T,S)
     @inbounds suma = convert(R, a[end])
@@ -28,24 +31,29 @@ function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, dx::S)
     end
     suma
 end
+
 evaluate{T<:Number}(a::Taylor1{T}) = a[1]
 
 doc"""
     evaluate(x, δt)
 
 Evaluates each element of `x::Array{Taylor1{T},1}`, representing
-the dependent variables of an ODE, at *time* δt.
+the dependent variables of an ODE, at *time* δt. Note that an array `x` of
+`Taylor1` polynomials may also be evaluated by calling it as a function;
+that is, the syntax `x(δt)` is equivalent to `evaluate(x, δt)`, and `x()`
+is equivalent to `evaluate(x)`.
 """
 function evaluate{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S)
     R = promote_type(T,S)
     return evaluate(convert(Array{Taylor1{R},1},x), convert(R,δt))
 end
+
 function evaluate{T<:Number}(x::Array{Taylor1{T},1}, δt::T)
     xnew = Array{T}( length(x) )
     evaluate!(x, δt, xnew)
     return xnew
 end
-evaluate{T<:Number}(a::Array{Taylor1{T},1}) = evaluate(a, zero(T))
+evaluate{T<:Number}(a::Array{Taylor1{T},1}) = evaluate.(a)
 
 doc"""
     evaluate!(x, δt, x0)
@@ -62,6 +70,7 @@ function evaluate!{T<:Number}(x::Array{Taylor1{T},1}, δt::T, x0::Union{Array{T,
     end
     nothing
 end
+
 function evaluate!{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S, x0::Union{Array{T,1},SubArray{T,1}})
     @assert length(x) == length(x0)
     @inbounds for i in eachindex(x)
@@ -70,11 +79,11 @@ function evaluate!{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S, x0::Uni
     nothing
 end
 
-
 """
     evaluate(a, x)
 
 Substitute `x::Taylor1` as independent variable in a `a::Taylor1` polynomial.
+Note that the syntax `a(x)` is equivalent to `evaluate(a, x)`.
 """
 evaluate{T<:Number,S<:Number}(a::Taylor1{T}, x::Taylor1{S}) =
     evaluate(promote(a,x)...)
@@ -89,7 +98,25 @@ function evaluate{T<:Number}(a::Taylor1{T}, x::Taylor1{T})
     suma
 end
 
+function evaluate{T<:Number}(a::Taylor1{Taylor1{T}}, x::Taylor1{T})
+    @inbounds suma = a[end]
+    @inbounds for k = a.order:-1:1
+        suma = suma*x + a[k]
+    end
+    suma
+end
+
+evaluate{T<:Number,S<:Number}(p::Taylor1{T},x::Array{S}) = evaluate.(p,x)
+
+#function-like behavior for Taylor1
+(p::Taylor1)(x) = evaluate(p, x)
+
+(p::Taylor1)() = evaluate(p)
+
+#function-like behavior for Array{Taylor1,1}
+(p::Array{Taylor1{T},1}){T<:Number}(x) = evaluate(p, x)
 
+(p::Array{Taylor1{T},1}){T<:Number}() = evaluate.(p)
 
 ## Evaluation of multivariable
 function evaluate!{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1},
@@ -100,6 +127,16 @@ function evaluate!{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1},
     end
     nothing
 end
+
+function evaluate!{T<:NumberNotSeriesN}(x::Array{TaylorN{T},1}, δx::Array{Taylor1{T},1},
+        x0::Array{Taylor1{T},1})
+    @assert length(x) == length(x0)
+    @inbounds for i in eachindex(x)
+        x0[i] = evaluate( x[i], δx )
+    end
+    nothing
+end
+
 function evaluate!{T<:Number}(x::Array{Taylor1{TaylorN{T}},1}, δt::T,
         x0::Array{TaylorN{T},1})
     @assert length(x) == length(x0)
@@ -110,9 +147,11 @@ function evaluate!{T<:Number}(x::Array{Taylor1{TaylorN{T}},1}, δt::T,
 end
 
 """
-    evaluate(a, vals)
+    evaluate(a, [vals])
 
-Evaluate a `HomogeneousPolynomial` polynomial at `vals`.
+Evaluate a `HomogeneousPolynomial` polynomial at `vals`. If `vals` is ommitted,
+it's evaluated at zero. Note that the syntax `a(vals)` is equivalent to
+`evaluate(a, vals)`; and `a()` is equivalent to `evaluate(a)`.
 """
 function evaluate{T<:Number,S<:NumberNotSeriesN}(a::HomogeneousPolynomial{T},
         vals::Array{S,1} )
@@ -133,14 +172,21 @@ function evaluate{T<:Number,S<:NumberNotSeriesN}(a::HomogeneousPolynomial{T},
 
     return suma
 end
+
 evaluate(a::HomogeneousPolynomial) = zero(a[1])
 
+#function-like behavior for HomogeneousPolynomial
+(p::HomogeneousPolynomial)(x) = evaluate(p, x)
+
+(p::HomogeneousPolynomial)() = evaluate(p)
 
 """
     evaluate(a, [vals])
 
 Evaluate the `TaylorN` polynomial `a` at `vals`.
 If `vals` is ommitted, it's evaluated at zero.
+Note that the syntax `a(vals)` is equivalent to `evaluate(a, vals)`; and `a()`
+is equivalent to `evaluate(a)`.
 """
 function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
         vals::Array{S,1})
@@ -165,6 +211,7 @@ function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
 
     return sum( sort!(suma, by=abs2) )
 end
+
 function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
         vals::Array{Taylor1{S},1})
     @assert length(vals) == get_numvars()
@@ -190,6 +237,30 @@ function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
     return suma
 end
 
+function evaluate{T<:NumberNotSeries}(a::TaylorN{Taylor1{T}},
+        vals::Array{Taylor1{T},1})
+    @assert length(vals) == get_numvars()
+
+    num_vars = get_numvars()
+    ct = coeff_table
+    a_length = length(a)
+    ord = maximum( get_order.(vals) )
+    suma = Taylor1(zeros(T, ord))
+
+    for homPol in 1:length(a)
+        sun = zero(Taylor1{T})
+        for (i,a_coeff) in enumerate(a.coeffs[homPol].coeffs)
+            tmp = vals[1]^(ct[homPol][i][1])
+            for n in 2:num_vars
+                tmp *= vals[n]^(ct[homPol][i][n])
+            end
+            suma += a_coeff * tmp
+        end
+    end
+
+    return suma
+end
+
 evaluate{T<:Number}(a::TaylorN{T}) = a[1][1]
 
 function evaluate{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1})
@@ -197,3 +268,22 @@ function evaluate{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1})
     evaluate!( x, δx, x0 )
     return x0
 end
+
+function evaluate{T<:NumberNotSeriesN}(x::Array{TaylorN{T},1}, δx::Array{Taylor1{T},1})
+    x0 = Array{Taylor1{T}}( length(x) )
+    evaluate!( x, δx, x0 )
+    return x0
+end
+
+evaluate{T<:Number}(x::Array{TaylorN{T},1}) = evaluate.(x)
+
+#function-like behavior for TaylorN
+(p::TaylorN)(x) = evaluate(p, x)
+
+(p::TaylorN)() = evaluate(p)
+
+#function-like behavior for Array{TaylorN,1}
+(p::Array{TaylorN{T},1}){T<:Number}(x) = evaluate(p, x)
+
+(p::Array{TaylorN{T},1}){T<:Number}() = evaluate(p)
+
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index 3ad22790..dce5ee8b 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -131,6 +131,12 @@ using Base.Test
     @test (xH/complex(0,BigInt(3)))' ==
         im*HomogeneousPolynomial([BigInt(1),0])/3
     @test evaluate(xH) == zero(eltype(xH))
+    @test xH() == zero(eltype(xH))
+    @test xH([1,1]) == evaluate(xH, [1,1])
+    hp = -5.4xH+6.89yH
+    @test hp([1,1]) == evaluate(hp, [1,1])
+    vr = rand(2)
+    @test hp(vr) == evaluate(hp, vr)
 
     @test derivative(2xT*yT^2,1) == 2yT^2
     @test xT*xT^3 == xT^4
@@ -341,4 +347,34 @@ using Base.Test
     @test_throws AssertionError cos(x)/sin(y)
     @test_throws BoundsError xH[20]
     @test_throws BoundsError xT[20]
+
+    #test function-like behavior for TaylorN
+    @test exy() == 1
+    @test exy([0.1im,0.01im]) == exp(0.11im)
+    @test isapprox(exy([1,1]), e^2)
+    @test sin(asin(xT+yT))([1.0,0.5]) == 1.5
+    @test asin(sin(xT+yT))([1.0,0.5]) == 1.5
+    @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(rand(2)) == 1
+    @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(zeros(2)) == 1
+    dx = set_variables("x", numvars=4, order=10)
+    P = sin.(dx)
+    v = [1.0,2,3,4]
+    for i in 1:4
+        @test P[i](v) == evaluate(P[i], v)
+    end
+    @test P.(fill(v, 4)) == fill(P(v), 4)
+    F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))]
+    Q = F(v+dx)
+    @test Q.( fill(v, 4) ) == fill(Q(v), 4)
+    vr = map(x->rand(4), 1:4)
+    @test Q.(vr) == map(x->Q(x), vr)
+    for i in 1:4
+        @test P[i]() == evaluate(P[i])
+        @test Q[i]() == evaluate(Q[i])
+    end
+    @test P() == evaluate.(P)
+    @test P() == evaluate(P)
+    @test Q() == evaluate.(Q)
+    @test Q() == evaluate(Q)
+    @test Q[1:3]() == evaluate(Q[1:3])
 end
diff --git a/test/mixtures.jl b/test/mixtures.jl
index 6ce96a9e..ff33177d 100644
--- a/test/mixtures.jl
+++ b/test/mixtures.jl
@@ -94,6 +94,7 @@ using Base.Test
     @test convert(Array{TaylorN{Taylor1{Float64}},1}, [t1N, t1N]) == [tN1, tN1]
     @test convert(Array{TaylorN{Taylor1{Float64}},2}, [t1N t1N]) == [tN1 tN1]
 
+    @test evaluate(t1N, 0.0) == TaylorN(xH, 2)
     @test string(evaluate(t1N, 0.0)) == " 1.0 x₁ + 𝒪(‖x‖³)"
     @test string(evaluate(t1N^2, 1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)"
     v = zeros(TaylorN{Float64},2)
@@ -133,4 +134,63 @@ using Base.Test
     @test xx*δx + Taylor1(typeof(δx),5) == δx + δx^2 + Taylor1(typeof(δx),5)
     @test xx/(1+δx) == one(xx)
     @test typeof(xx+δx) == Taylor1{TaylorN{Float64}}
+
+    #testing evaluate and function-like behavior of Taylor1, TaylorN for mixtures:
+    t = Taylor1(25)
+    p = cos(t)
+    q = sin(t)
+    a = [p,q]
+    dx = set_variables("x", numvars=4, order=10)
+    P = sin.(dx)
+    v = [1.0,2,3,4]
+    F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))]
+    Q = F(v+dx)
+    diff_evals = cos(sin(dx[1]))-p(P[1])
+    @test norm( norm.(map(x->x.coeffs, diff_evals.coeffs),Inf) , Inf) < 1e-15
+    #evaluate a Taylor1 at a TaylorN
+    @test p(P) == evaluate(p, P)
+    @test q(Q) == evaluate(q, Q)
+    #evaluate an array of Taylor1s at a TaylorN
+    aT1 = [p,q,p^2,log(1+q)] #an array of Taylor1s
+    @test aT1(Q[4]) == evaluate(aT1, Q[4])
+    @test (aT1.^2)(Q[3]) == evaluate(aT1.^2, Q[3])
+    #evaluate a TaylorN at an array of Taylor1s
+    @test P[1](aT1) == evaluate(P[1], aT1)
+    @test Q[2](aT1) == evaluate(Q[2], aT1)
+    #evaluate an array of TaylorN{Float64} at an array of Taylor1{Float64}
+    @test P(aT1) == evaluate(P, aT1)
+    @test Q(aT1) == evaluate(Q, aT1)
+    #test evaluation of an Array{TaylorN{Taylor1}} at an Array{Taylor1}
+    aH1 = [
+        HomogeneousPolynomial([Taylor1(rand(2))]), 
+        HomogeneousPolynomial([Taylor1(rand(2)),Taylor1(rand(2)),Taylor1(rand(2)),Taylor1(rand(2))])
+        ]
+    bH1 = [
+        HomogeneousPolynomial([Taylor1(rand(2))]), 
+        HomogeneousPolynomial([Taylor1(rand(2)),Taylor1(rand(2)),Taylor1(rand(2)),Taylor1(rand(2))])
+        ]
+    aTN1 = TaylorN(aH1); bTN1 = TaylorN(bH1)
+    x = [aTN1, bTN1]
+    δx = [Taylor1(rand(3)) for i in 1:4]
+    @test typeof(x) == Array{TaylorN{Taylor1{Float64}},1}
+    @test typeof(δx) == Array{Taylor1{Float64},1}
+    x0 = Array{Taylor1{Float64}}(length(x))
+    eval_x_δx = evaluate(x,δx)
+    @test x(δx) == evaluate(x,δx)
+    evaluate!(x,δx,x0)
+    @test x0 == eval_x_δx
+    @test typeof(evaluate(x[1],δx)) == Taylor1{Float64}
+    @test x() == map(y->y[1][1], x)
+    for i in eachindex(x)
+        @test evaluate(x[i],δx) == eval_x_δx[i]
+        @test x[i](δx) == eval_x_δx[i]
+    end
+    p11 = Taylor1([sin(t),cos(t)])
+    @which evaluate(p11,t)
+    @test evaluate(p11,t) == sin(t)+t*cos(t)
+    @test p11(t) == sin(t)+t*cos(t)
+    a11 = Taylor1([t,t^2,exp(-t),sin(t),cos(t)])
+    b11 = t+t*(t^2)+(t^2)*(exp(-t))+(t^3)*sin(t)+(t^4)*cos(t)
+    diff_a11b11 = a11(t)-b11
+    @test norm(diff_a11b11.coeffs,Inf) < 1E-19
 end
diff --git a/test/onevariable.jl b/test/onevariable.jl
index e102b8bf..ceb6fc5e 100644
--- a/test/onevariable.jl
+++ b/test/onevariable.jl
@@ -164,6 +164,48 @@ using Base.Test
     @test evaluate(exp(Taylor1([0,1],17)),1.0) == 1.0*e
     @test evaluate(exp(Taylor1([0,1],1))) == 1.0
     @test evaluate(exp(t),t^2) == exp(t^2)
+    #Test function-like behavior for Taylor1s
+    t17 = Taylor1([0,1],17)
+    myexpfun = exp(t17)
+    @test myexpfun(1.0) == 1.0*e
+    @test myexpfun() == 1.0
+    @test myexpfun(t17^2) == exp(t17^2)
+    @test exp(t17^2)(t17) == exp(t17^2)
+    p = cos(t17)
+    q = sin(t17)
+    @test cos(-im*t)(1)+im*sin(-im*t)(1) == exp(-im*t)(im)
+    @test p(-im*t17)(1)+im*q(-im*t17)(1) == exp(-im*t17)(im)
+    cossin1 = x->p(q(x))
+    @test evaluate(p, evaluate(q, pi/4)) == cossin1(pi/4)
+    cossin2 = p(q)
+    @test evaluate(evaluate(p,q), pi/4) == cossin2(pi/4)
+    @test evaluate(p, q) == cossin2
+    @test p(q)() == evaluate(evaluate(p, q))
+    @test evaluate(p, q) == p(q)
+    @test evaluate(q, p) == q(p)
+    cs = x->cos(sin(x))
+    csdiff = (cs(t17)-cossin2(t17)).(-2:0.1:2)
+    @test norm(csdiff, 1) < 5e-15
+    a = [p,q]
+    @test a(0.1) == evaluate.([p,q],0.1)
+    @test a.(0.1) == a(0.1)
+    @test a.() == evaluate.([p, q])
+    @test a.() == [p(), q()]
+    @test a.() == a()
+    vr = rand(2)
+    @test p.(vr) == evaluate.(p, vr)
+    Mr = rand(3,3,3)
+    @test p.(Mr) == evaluate.(p, Mr)
+    mytaylor1 = Taylor1(rand(20))
+    vr = rand(5)
+    @test p(vr) == p.(vr)
+    @test p(vr) == evaluate.(p,vr)
+    @test p(Mr) == p.(Mr)
+    @test p(Mr) == evaluate.(p,Mr)
+    taylor_a = Taylor1(Int64,10)
+    taylor_x = exp(Taylor1(Float64,13))
+    @which evaluate(taylor_x, taylor_a)
+    @test taylor_x(taylor_a) == evaluate(taylor_x, taylor_a)
 
     @test sin(asin(tsquare)) == tsquare
     @test tan(atan(tsquare)) == tsquare