diff --git a/include/dacex/Taylor.h b/include/dacex/Taylor.h
index d4e666d..0e301fe 100644
--- a/include/dacex/Taylor.h
+++ b/include/dacex/Taylor.h
@@ -198,6 +198,47 @@ class Polynomial<Taylor> : public PolynomialBase {
      */
     friend Polynomial<Taylor> operator/(const double c, const Polynomial<Taylor> &p);
 
+    /** @brief Evaluate at polynomial-vector arguments.
+     * @param args Algebraic-vector argument values.
+     * @return Evaluated Taylor polynomial.
+     */
+    Polynomial<Taylor> eval(const DACE::AlgebraicVector<Polynomial<Taylor>>& args) const;
+    /** @brief Evaluate at numeric algebraic-vector arguments.
+     * @param args Numeric argument values.
+     * @return Evaluated scalar value.
+     */
+    double eval(const DACE::AlgebraicVector<double>& args) const;
+    /** @brief Evaluate at numeric vector arguments.
+     * @param args Numeric argument values.
+     * @return Evaluated scalar value.
+     */
+    double eval(const std::vector<double>& args) const;
+    /** @brief Evaluate with generic vector arguments.
+     * @tparam T Evaluation scalar type.
+     * @param args Input argument vector.
+     * @return Evaluated value in type T.
+     */
+    template<class T> T eval(const std::vector<T> &args) const {
+        return m_da.eval(args);
+    }
+    /** @brief Evaluate with generic array arguments.
+     * @tparam T Evaluation scalar type.
+     * @param args Input argument array.
+     * @param length Number of array elements.
+     * @return Evaluated value in type T.
+     */
+    template<class T> T eval(const T args[], const unsigned int length) const {
+        return m_da.eval(args, length);
+    }
+    /** @brief Evaluate a single-variable polynomial.
+     * @tparam T Evaluation scalar type.
+     * @param arg Input scalar argument.
+     * @return Evaluated value in type T.
+     */
+    template<class T> T evalScalar(const T &arg) const {
+        return m_da.evalScalar(arg);
+    }
+
     /** @brief Multiply matching monomials coefficient-wise.
      * @param p Right-hand side polynomial.
      * @return Product polynomial.
diff --git a/src/Taylor.cpp b/src/Taylor.cpp
index a619ecc..8948aca 100644
--- a/src/Taylor.cpp
+++ b/src/Taylor.cpp
@@ -31,6 +31,26 @@ DACE::AlgebraicVector<Polynomial<Taylor>> Polynomial<Taylor>::gradient() const {
     return result;
 }
 
+Polynomial<Taylor> Polynomial<Taylor>::eval(const DACE::AlgebraicVector<Polynomial<Taylor>>& args) const {
+    std::vector<DACE::DA> da_args(args.size());
+    for (size_t i = 0; i < args.size(); ++i) {
+        da_args[i] = *args[i].getDA();
+    }
+    return Polynomial<Taylor>(m_da.eval(da_args));
+}
+
+double Polynomial<Taylor>::eval(const DACE::AlgebraicVector<double>& args) const {
+    std::vector<double> std_args(args.size());
+    for (size_t i = 0; i < args.size(); ++i) {
+        std_args[i] = args[i];
+    }
+    return m_da.eval(std_args);
+}
+
+double Polynomial<Taylor>::eval(const std::vector<double>& args) const {
+    return m_da.eval(args);
+}
+
 
 
 Polynomial<Taylor>& Polynomial<Taylor>::operator=(Polynomial<Taylor> &&p) noexcept {