diff --git a/stan/math/prim/err/check_consistent_sizes.hpp b/stan/math/prim/err/check_consistent_sizes.hpp
index ec1cee072e6..d01636e33b0 100644
--- a/stan/math/prim/err/check_consistent_sizes.hpp
+++ b/stan/math/prim/err/check_consistent_sizes.hpp
@@ -1,8 +1,10 @@
 #ifndef STAN_MATH_PRIM_ERR_CHECK_CONSISTENT_SIZES_HPP
 #define STAN_MATH_PRIM_ERR_CHECK_CONSISTENT_SIZES_HPP
 
+#include <stan/math/prim/err/check_matching_sizes.hpp>
 #include <stan/math/prim/err/invalid_argument.hpp>
 #include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/meta/is_container.hpp>
 #include <stan/math/prim/meta/require_generics.hpp>
 #include <algorithm>
 #include <sstream>
@@ -66,6 +68,34 @@ inline void check_consistent_sizes(const char* function, const char* name1,
   }
 }
 
+/**
+ * Check that the unnamed container inputs have the same size. Inputs are
+ * labeled `arg1`, `arg2`, and so on in error messages.
+ *
+ * @tparam T type of the first input
+ * @tparam Types types of the remaining inputs
+ * @param function function name (for error messages)
+ * @param x first input
+ * @param xs remaining inputs
+ * @throw `invalid_argument` if sizes are inconsistent
+ */
+template <typename T, typename... Types,
+          require_all_container_t<T, Types...>* = nullptr>
+inline void check_consistent_sizes(const char* function, T&& x, Types&&... xs) {
+  std::size_t arg_idx = 2;
+  (
+      [&](const auto& y) {
+        if (x.size() != y.size()) {
+          [&]() STAN_COLD_PATH {
+            const std::string name = "arg" + std::to_string(arg_idx);
+            check_matching_sizes(function, "arg1", x, name.c_str(), y);
+          }();
+        }
+        ++arg_idx;
+      }(xs),
+      ...);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/err/check_matching_dims.hpp b/stan/math/prim/err/check_matching_dims.hpp
index 5308545834f..7856b9a7b97 100644
--- a/stan/math/prim/err/check_matching_dims.hpp
+++ b/stan/math/prim/err/check_matching_dims.hpp
@@ -155,6 +155,35 @@ inline void check_matching_dims(const char* function, const char* name1,
   check_matching_dims(function, name1, y1, name2, y2);
 }
 
+/**
+ * Check that the unnamed Eigen inputs have the same dimensions. Inputs are
+ * labeled `arg1`, `arg2`, and so on in error messages.
+ *
+ * @tparam T type of the first input
+ * @tparam Types types of the remaining inputs
+ * @param function function name (for error messages)
+ * @param x first input
+ * @param xs remaining inputs
+ * @throw `invalid_argument` if dimensions do not match
+ */
+template <typename T, typename... Types,
+          require_all_eigen_t<T, Types...>* = nullptr>
+inline void check_matching_dims(const char* function, const T& x,
+                                const Types&... xs) {
+  std::size_t arg_idx = 2;
+  (
+      [&](const auto& y) {
+        if (x.rows() != y.rows() || x.cols() != y.cols()) {
+          [&]() STAN_COLD_PATH {
+            const std::string name = "arg" + std::to_string(arg_idx);
+            check_matching_dims(function, "arg1", x, name.c_str(), y);
+          }();
+        }
+        ++arg_idx;
+      }(xs),
+      ...);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/functor.hpp b/stan/math/prim/functor.hpp
index 2f9de59c706..22a3934f99d 100644
--- a/stan/math/prim/functor.hpp
+++ b/stan/math/prim/functor.hpp
@@ -23,6 +23,7 @@
 #include <stan/math/prim/functor/ode_ckrk.hpp>
 #include <stan/math/prim/functor/ode_rk45.hpp>
 #include <stan/math/prim/functor/ode_store_sensitivities.hpp>
+#include <stan/math/prim/functor/map.hpp>
 #include <stan/math/prim/functor/map_if.hpp>
 #include <stan/math/prim/functor/map_rect.hpp>
 #include <stan/math/prim/functor/map_rect_combine.hpp>
diff --git a/stan/math/prim/functor/map.hpp b/stan/math/prim/functor/map.hpp
new file mode 100644
index 00000000000..0eb5740b733
--- /dev/null
+++ b/stan/math/prim/functor/map.hpp
@@ -0,0 +1,384 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_MAP_HPP
+#define STAN_MATH_PRIM_FUNCTOR_MAP_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err/check_consistent_sizes.hpp>
+#include <stan/math/prim/err/check_matching_dims.hpp>
+#include <stan/math/prim/err/check_size_match.hpp>
+#include <stan/math/prim/functor/apply.hpp>
+
+#include <cstddef>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return a std::vector whose i-th element is `f(x[i], args...)`.
+ * Additional arguments are shared across calls.  If `x` is empty, the result
+ * is empty.  The element type of the result is deduced from `f`.
+ *
+ * @tparam F type of functor
+ * @tparam T type of input std::vector
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each element
+ * @param[in] x input std::vector
+ * @param[in] args shared arguments passed to each call
+ * @return std::vector of results
+ */
+template <typename F, typename T, typename... Args,
+          require_std_vector_t<T>* = nullptr>
+inline auto map(F&& f, T&& x, Args&&... args) {
+  using T_return = std::decay_t<decltype(f(x[0], args...))>;
+  std::vector<T_return> result;
+  result.reserve(x.size());
+
+  for (auto&& xi : x) {
+    result.push_back(f(xi, args...));
+  }
+
+  return result;
+}
+
+/**
+ * Return a std::vector whose i-th element is
+ * `f(get<0>(xs)[i], get<1>(xs)[i], ..., args...)`.
+ * The tuple `xs` must contain at least two equal-sized std::vector inputs.
+ * Additional trailing arguments are shared across calls.  If the inputs are
+ * empty, the result is empty.  The element type of the result is deduced from
+ * `f`.
+ *
+ * @tparam F type of functor
+ * @tparam Tuple type of tuple of std::vector inputs
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each tuple of elements
+ * @param[in] xs tuple of std::vector inputs
+ * @param[in] args shared arguments passed to each call
+ * @return std::vector of results
+ * @throw std::invalid_argument if the inputs have different sizes
+ */
+template <typename F, typename Tuple, typename... Args,
+          require_tuple_t<Tuple>* = nullptr>
+inline auto mapN(F&& f, Tuple&& xs, Args&&... args) {
+  static_assert(tuple_size_v<Tuple> >= 2,
+                "mapN requires a tuple of at least two std::vector inputs.");
+  static constexpr const char* function = "mapN";
+
+  return apply(
+      [&](auto&&... xs) {
+        static_assert((is_std_vector<std::decay_t<decltype(xs)>>::value && ...),
+                      "mapN tuple elements must be std::vector.");
+        check_consistent_sizes(function, xs...);
+
+        const std::size_t n = std::get<0>(std::forward_as_tuple(xs...)).size();
+        using T_return = std::decay_t<decltype(f((xs[0])..., args...))>;
+        std::vector<T_return> result;
+        result.reserve(n);
+
+        for (std::size_t i = 0; i < n; ++i) {
+          result.push_back(f((xs[i])..., args...));
+        }
+        return result;
+      },
+      std::forward<Tuple>(xs));
+}
+
+/**
+ * Return a std::vector whose i-th element is `f(x1[i], x2[i], ...)`.
+ * Convenience overload equivalent to
+ * `mapN(f, std::forward_as_tuple(x1, x2, xs...))`.
+ *
+ * @tparam F type of functor
+ * @tparam T1 type of first std::vector
+ * @tparam T2 type of second std::vector
+ * @tparam Ts types of additional std::vector inputs
+ * @param[in] f functor applied to each tuple of elements
+ * @param[in] x1 first std::vector
+ * @param[in] x2 second std::vector
+ * @param[in] xs additional std::vector inputs
+ * @return std::vector of results
+ * @throw std::invalid_argument if the inputs have different sizes
+ */
+template <typename F, typename T1, typename T2, typename... Ts,
+          require_all_std_vector_t<T1, T2, Ts...>* = nullptr>
+inline auto mapN(F&& f, T1&& x1, T2&& x2, Ts&&... xs) {
+  return mapN(std::forward<F>(f), std::forward_as_tuple(x1, x2, xs...));
+}
+
+/**
+ * Return a matrix whose i-th row is `f(m.row(i), args...)`.
+ * Additional arguments are shared across calls.  If `m` has zero rows, or if
+ * the first returned row is empty, returns a 0x0 matrix.  All returned rows
+ * must have the same number of columns.
+ *
+ * @tparam F type of functor
+ * @tparam T type of Eigen matrix
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each row
+ * @param[in] m input matrix
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if returned rows have inconsistent sizes
+ */
+template <typename F, typename T, typename... Args,
+          require_eigen_matrix_dynamic_t<T>* = nullptr>
+inline auto row_map(F&& f, T&& m, Args&&... args) {
+  static constexpr const char* function = "row_map";
+  decltype(auto) m_ref = to_ref(std::forward<T>(m));
+  using result_row_t = plain_type_t<decltype(f(m_ref.row(0), args...))>;
+  using T_return = scalar_type_t<result_row_t>;
+  using matrix_t = Eigen::Matrix<T_return, Eigen::Dynamic, Eigen::Dynamic>;
+
+  const Eigen::Index n_rows = m_ref.rows();
+  if (n_rows == 0) {
+    return matrix_t(0, 0);
+  }
+
+  result_row_t first_row = f(m_ref.row(0), args...);
+  if (first_row.size() == 0) {
+    return matrix_t(0, 0);
+  }
+  matrix_t result(n_rows, first_row.cols());
+  result.row(0) = first_row;
+  for (Eigen::Index i = 1; i < n_rows; ++i) {
+    result_row_t row = f(m_ref.row(i), args...);
+    check_size_match(function, "columns of result", result.cols(),
+                     "columns of returned row", row.cols());
+    result.row(i) = row;
+  }
+  return result;
+}
+
+/**
+ * Return a matrix whose j-th column is `f(m.col(j), args...)`.
+ * Additional arguments are shared across calls.  If `m` has zero columns, or
+ * if the first returned column is empty, returns a 0x0 matrix.  All returned
+ * columns must have the same number of rows.
+ *
+ * @tparam F type of functor
+ * @tparam T type of Eigen matrix
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each column
+ * @param[in] m input matrix
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if returned columns have inconsistent sizes
+ */
+template <typename F, typename T, typename... Args,
+          require_eigen_matrix_dynamic_t<T>* = nullptr>
+inline auto col_map(F&& f, T&& m, Args&&... args) {
+  static constexpr const char* function = "col_map";
+  decltype(auto) m_ref = to_ref(std::forward<T>(m));
+  using result_col_t = plain_type_t<decltype(f(m_ref.col(0), args...))>;
+  using T_return = scalar_type_t<result_col_t>;
+  using matrix_t = Eigen::Matrix<T_return, Eigen::Dynamic, Eigen::Dynamic>;
+
+  const Eigen::Index n_cols = m_ref.cols();
+  if (n_cols == 0) {
+    return matrix_t(0, 0);
+  }
+
+  result_col_t first_col = f(m_ref.col(0), args...);
+  if (first_col.size() == 0) {
+    return matrix_t(0, 0);
+  }
+  matrix_t result(first_col.rows(), n_cols);
+  result.col(0) = first_col;
+  for (Eigen::Index j = 1; j < n_cols; ++j) {
+    result_col_t col = f(m_ref.col(j), args...);
+    check_size_match(function, "rows of result", result.rows(),
+                     "rows of returned column", col.rows());
+    result.col(j) = col;
+  }
+  return result;
+}
+
+/**
+ * Return a matrix whose i-th row is
+ * `f(get<0>(ms).row(i), get<1>(ms).row(i), ..., args...)`.
+ * The tuple `ms` must contain at least two Eigen matrices with identical
+ * dimensions.  Additional trailing arguments are shared across calls.  If the
+ * inputs have zero rows, or if the first returned row is empty, returns a 0x0
+ * matrix.  All returned rows must have the same number of columns.
+ *
+ * @tparam F type of functor
+ * @tparam Tuple type of tuple of Eigen matrix inputs
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each tuple of rows
+ * @param[in] ms tuple of Eigen matrix inputs
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if the inputs have different dimensions or
+ * returned rows have inconsistent sizes
+ */
+template <typename F, typename Tuple, typename... Args,
+          require_tuple_t<Tuple>* = nullptr>
+inline auto row_mapN(F&& f, Tuple&& ms, Args&&... args) {
+  static_assert(tuple_size_v<Tuple> >= 2,
+                "row_mapN requires a tuple of at least two Eigen matrix "
+                "inputs.");
+  static constexpr const char* function = "row_mapN";
+
+  return apply(
+      [&](auto&&... mats) {
+        static_assert(
+            (is_eigen_matrix_dynamic<std::decay_t<decltype(mats)>>::value
+             && ...),
+            "row_mapN tuple elements must be Eigen matrices.");
+        check_matching_dims(function, mats...);
+
+        auto m_refs = std::tuple{to_ref(std::forward<decltype(mats)>(mats))...};
+        const Eigen::Index n_rows = std::get<0>(m_refs).rows();
+        using result_row_t = plain_type_t<decltype(f(
+            (std::declval<ref_type_t<decltype(mats)>>().row(0))..., args...))>;
+        using T_return = scalar_type_t<result_row_t>;
+        using matrix_t
+            = Eigen::Matrix<T_return, Eigen::Dynamic, Eigen::Dynamic>;
+
+        if (n_rows == 0) {
+          return matrix_t(0, 0);
+        }
+
+        return apply(
+            [&](auto&&... mrs) {
+              result_row_t first_row = f((mrs.row(0))..., args...);
+              if (first_row.size() == 0) {
+                return matrix_t(0, 0);
+              }
+              matrix_t result(n_rows, first_row.cols());
+              result.row(0) = first_row;
+              for (Eigen::Index i = 1; i < n_rows; ++i) {
+                result_row_t row = f((mrs.row(i))..., args...);
+                check_size_match(function, "columns of result", result.cols(),
+                                 "columns of returned row", row.cols());
+                result.row(i) = row;
+              }
+              return result;
+            },
+            m_refs);
+      },
+      std::forward<Tuple>(ms));
+}
+
+/**
+ * Return a matrix whose i-th row is `f(m1.row(i), m2.row(i), ...)`.
+ * Convenience overload equivalent to
+ * `row_mapN(f, std::forward_as_tuple(m1, m2, ms...))`.
+ *
+ * @tparam F type of functor
+ * @tparam T1 type of first Eigen matrix
+ * @tparam T2 type of second Eigen matrix
+ * @tparam Ts types of additional Eigen matrix inputs
+ * @param[in] f functor applied to each tuple of rows
+ * @param[in] m1 first Eigen matrix
+ * @param[in] m2 second Eigen matrix
+ * @param[in] ms additional Eigen matrix inputs
+ * @return matrix of results
+ * @throw std::invalid_argument if the inputs have different dimensions or
+ * returned rows have inconsistent sizes
+ */
+template <typename F, typename T1, typename T2, typename... Ts,
+          require_all_eigen_matrix_dynamic_t<T1, T2, Ts...>* = nullptr>
+inline auto row_mapN(F&& f, T1&& m1, T2&& m2, Ts&&... ms) {
+  return row_mapN(std::forward<F>(f), std::forward_as_tuple(m1, m2, ms...));
+}
+
+/**
+ * Return a matrix whose j-th column is
+ * `f(get<0>(ms).col(j), get<1>(ms).col(j), ..., args...)`.
+ * The tuple `ms` must contain at least two Eigen matrices with identical
+ * dimensions.  Additional trailing arguments are shared across calls.  If the
+ * inputs have zero columns, or if the first returned column is empty, returns
+ * a 0x0 matrix.  All returned columns must have the same number of rows.
+ *
+ * @tparam F type of functor
+ * @tparam Tuple type of tuple of Eigen matrix inputs
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each tuple of columns
+ * @param[in] ms tuple of Eigen matrix inputs
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if the inputs have different dimensions or
+ * returned columns have inconsistent sizes
+ */
+template <typename F, typename Tuple, typename... Args,
+          require_tuple_t<Tuple>* = nullptr>
+inline auto col_mapN(F&& f, Tuple&& ms, Args&&... args) {
+  static_assert(tuple_size_v<Tuple> >= 2,
+                "col_mapN requires a tuple of at least two Eigen matrix "
+                "inputs.");
+  static constexpr const char* function = "col_mapN";
+
+  return apply(
+      [&](auto&&... mats) {
+        static_assert(
+            (is_eigen_matrix_dynamic<std::decay_t<decltype(mats)>>::value
+             && ...),
+            "col_mapN tuple elements must be Eigen matrices.");
+        check_matching_dims(function, mats...);
+
+        auto m_refs = std::tuple{to_ref(std::forward<decltype(mats)>(mats))...};
+        const Eigen::Index n_cols = std::get<0>(m_refs).cols();
+        using result_col_t = plain_type_t<decltype(f(
+            (std::declval<ref_type_t<decltype(mats)>>().col(0))..., args...))>;
+        using T_return = scalar_type_t<result_col_t>;
+        using matrix_t
+            = Eigen::Matrix<T_return, Eigen::Dynamic, Eigen::Dynamic>;
+
+        if (n_cols == 0) {
+          return matrix_t(0, 0);
+        }
+
+        return apply(
+            [&](auto&&... mrs) {
+              result_col_t first_col = f((mrs.col(0))..., args...);
+              if (first_col.size() == 0) {
+                return matrix_t(0, 0);
+              }
+              matrix_t result(first_col.rows(), n_cols);
+              result.col(0) = first_col;
+              for (Eigen::Index j = 1; j < n_cols; ++j) {
+                result_col_t col = f((mrs.col(j))..., args...);
+                check_size_match(function, "rows of result", result.rows(),
+                                 "rows of returned column", col.rows());
+                result.col(j) = col;
+              }
+              return result;
+            },
+            m_refs);
+      },
+      std::forward<Tuple>(ms));
+}
+
+/**
+ * Return a matrix whose j-th column is `f(m1.col(j), m2.col(j), ...)`.
+ * Convenience overload equivalent to
+ * `col_mapN(f, std::forward_as_tuple(m1, m2, ms...))`.
+ *
+ * @tparam F type of functor
+ * @tparam T1 type of first Eigen matrix
+ * @tparam T2 type of second Eigen matrix
+ * @tparam Ts types of additional Eigen matrix inputs
+ * @param[in] f functor applied to each tuple of columns
+ * @param[in] m1 first Eigen matrix
+ * @param[in] m2 second Eigen matrix
+ * @param[in] ms additional Eigen matrix inputs
+ * @return matrix of results
+ * @throw std::invalid_argument if the inputs have different dimensions or
+ * returned columns have inconsistent sizes
+ */
+template <typename F, typename T1, typename T2, typename... Ts,
+          require_all_eigen_matrix_dynamic_t<T1, T2, Ts...>* = nullptr>
+inline auto col_mapN(F&& f, T1&& m1, T2&& m2, Ts&&... ms) {
+  return col_mapN(std::forward<F>(f), std::forward_as_tuple(m1, m2, ms...));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/prim/err/check_consistent_sizes_test.cpp b/test/unit/math/prim/err/check_consistent_sizes_test.cpp
new file mode 100644
index 00000000000..62fcc9eaa46
--- /dev/null
+++ b/test/unit/math/prim/err/check_consistent_sizes_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/prim/err/check_consistent_sizes.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+#include <vector>
+
+TEST(ErrorHandling, checkConsistentSizesUnnamedArguments) {
+  using stan::math::check_consistent_sizes;
+
+  std::vector<double> x1(3);
+  std::vector<int> x2(3);
+  std::vector<double> x3(2);
+
+  EXPECT_NO_THROW(check_consistent_sizes("check_consistent_sizes", x1, x2));
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_consistent_sizes("check_consistent_sizes", x1, x2, x3),
+      std::invalid_argument, "arg1", 1);
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_consistent_sizes("check_consistent_sizes", x1, x2, x3),
+      std::invalid_argument, "arg3", 1);
+  EXPECT_NO_THROW(check_consistent_sizes("check_consistent_sizes", x1));
+}
+
+TEST(ErrorHandling, checkConsistentSizesNamedArgumentsStillResolve) {
+  using stan::math::check_consistent_sizes;
+
+  std::vector<double> x1(3);
+  std::vector<double> x2(2);
+
+  EXPECT_THROW_MSG(
+      check_consistent_sizes("check_consistent_sizes", "x1", x1, "x2", x2),
+      std::invalid_argument, "x2");
+}
diff --git a/test/unit/math/prim/err/check_matching_dims_test.cpp b/test/unit/math/prim/err/check_matching_dims_test.cpp
index 4643355b52f..f0b62a70493 100644
--- a/test/unit/math/prim/err/check_matching_dims_test.cpp
+++ b/test/unit/math/prim/err/check_matching_dims_test.cpp
@@ -27,6 +27,23 @@ TEST(ErrorHandlingMatrix, checkMatchingDimsMatrix) {
       std::invalid_argument);
 }
 
+TEST(ErrorHandlingMatrix, checkMatchingDimsUnnamedArguments) {
+  using stan::math::check_matching_dims;
+
+  Eigen::MatrixXd x1(2, 3);
+  Eigen::MatrixXd x2(2, 3);
+  Eigen::MatrixXd x3(3, 2);
+
+  EXPECT_NO_THROW(check_matching_dims("checkMatchingDims", x1));
+  EXPECT_NO_THROW(check_matching_dims("checkMatchingDims", x1, x2));
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_matching_dims("checkMatchingDims", x1, x2, x3),
+      std::invalid_argument, "arg1", 1);
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_matching_dims("checkMatchingDims", x1, x2, x3),
+      std::invalid_argument, "arg3", 1);
+}
+
 TEST(ErrorHandlingMatrix, checkMatchingDimsArray) {
   std::vector<double> y(3);
   std::vector<double> x(3);
diff --git a/test/unit/math/prim/functor/map_test.cpp b/test/unit/math/prim/functor/map_test.cpp
new file mode 100644
index 00000000000..f4ddeee9eae
--- /dev/null
+++ b/test/unit/math/prim/functor/map_test.cpp
@@ -0,0 +1,531 @@
+#include <stan/math/prim.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <memory>
+#include <stdexcept>
+#include <tuple>
+#include <vector>
+
+TEST(MathFunctor, map_square) {
+  std::vector<double> x{1.0, 2.0, 3.0};
+  auto y = stan::math::map([](double v) { return v * v; }, x);
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 1.0);
+  EXPECT_FLOAT_EQ(y[1], 4.0);
+  EXPECT_FLOAT_EQ(y[2], 9.0);
+}
+
+TEST(MathFunctor, map_type_change) {
+  std::vector<int> x{1, 2, 3};
+  auto y = stan::math::map([](int v) { return v + 0.5; }, x);
+  EXPECT_FLOAT_EQ(y[0], 1.5);
+  EXPECT_FLOAT_EQ(y[1], 2.5);
+  EXPECT_FLOAT_EQ(y[2], 3.5);
+}
+
+TEST(MathFunctor, map_empty) {
+  std::vector<double> x;
+  auto y = stan::math::map([](double v) { return v; }, x);
+  EXPECT_EQ(y.size(), 0u);
+}
+
+TEST(MathFunctor, map_const_vector) {
+  const std::vector<double> x{1.0, 2.0, 3.0};
+  auto y = stan::math::map([](double v) { return 2.0 * v; }, x);
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 2.0);
+  EXPECT_FLOAT_EQ(y[1], 4.0);
+  EXPECT_FLOAT_EQ(y[2], 6.0);
+}
+
+TEST(MathFunctor, map_variadic_args) {
+  std::vector<double> x{1.0, 2.0, 3.0};
+  double offset = 10.0;
+  double scale = 2.0;
+  auto y = stan::math::map(
+      [](double v, double mu, double sigma) { return mu + sigma * v; }, x,
+      offset, scale);
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 12.0);
+  EXPECT_FLOAT_EQ(y[1], 14.0);
+  EXPECT_FLOAT_EQ(y[2], 16.0);
+}
+
+TEST(MathFunctor, map_shared_arg_reused) {
+  std::vector<double> x{1.0, 2.0, 3.0};
+  auto offset = std::make_unique<double>(10.0);
+  auto y = stan::math::map(
+      [](double v, const std::unique_ptr<double>& mu) { return v + *mu; }, x,
+      std::move(offset));
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 11.0);
+  EXPECT_FLOAT_EQ(y[1], 12.0);
+  EXPECT_FLOAT_EQ(y[2], 13.0);
+}
+
+TEST(MathFunctor, mapN_add) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0, 20.0};
+  auto y = stan::math::mapN([](double u, double v) { return u + v; }, a, b);
+  ASSERT_EQ(y.size(), 2u);
+  EXPECT_FLOAT_EQ(y[0], 11.0);
+  EXPECT_FLOAT_EQ(y[1], 22.0);
+}
+
+TEST(MathFunctor, mapN_mixed_scalar_types) {
+  std::vector<int> a{1, 2, 3};
+  std::vector<double> b{0.5, 1.5, 2.5};
+  auto y = stan::math::mapN([](int u, double v) { return u + v; }, a, b);
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 1.5);
+  EXPECT_FLOAT_EQ(y[1], 3.5);
+  EXPECT_FLOAT_EQ(y[2], 5.5);
+}
+
+TEST(MathFunctor, mapN_three_vectors) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0, 20.0};
+  std::vector<double> c{100.0, 200.0};
+  auto y = stan::math::mapN(
+      [](double u, double v, double w) { return u + v + w; }, a, b, c);
+  ASSERT_EQ(y.size(), 2u);
+  EXPECT_FLOAT_EQ(y[0], 111.0);
+  EXPECT_FLOAT_EQ(y[1], 222.0);
+}
+
+TEST(MathFunctor, mapN_shared_arg_reused) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0, 20.0};
+  auto offset = std::make_unique<double>(100.0);
+  auto y = stan::math::mapN(
+      [](double u, double v, const std::unique_ptr<double>& mu) {
+        return u + v + *mu;
+      },
+      std::forward_as_tuple(a, b), std::move(offset));
+  EXPECT_STD_VECTOR_FLOAT_EQ(y, std::vector<double>({111.0, 122.0}));
+}
+
+TEST(MathFunctor, mapN_empty) {
+  std::vector<double> a;
+  std::vector<double> b;
+  auto y = stan::math::mapN([](double u, double v) { return u + v; }, a, b);
+  EXPECT_EQ(y.size(), 0u);
+}
+
+TEST(MathFunctor, mapN_size_mismatch) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0};
+  EXPECT_THROW(stan::math::mapN([](double u, double v) { return u + v; }, a, b),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, row_map_standardize) {
+  Eigen::MatrixXd x(2, 3);
+  x << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0;
+  Eigen::RowVectorXd mu(3);
+  mu << 1.0, 2.0, 3.0;
+  Eigen::RowVectorXd sigma(3);
+  sigma << 1.0, 2.0, 3.0;
+
+  auto y = stan::math::row_map(
+      [](const Eigen::RowVectorXd& row, const Eigen::RowVectorXd& m,
+         const Eigen::RowVectorXd& s) {
+        return (row.array() - m.array()) / s.array();
+      },
+      x, mu, sigma);
+
+  Eigen::MatrixXd expected(2, 3);
+  expected << 0.0, 0.0, 0.0, 3.0, 1.5, 1.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_map_empty) {
+  Eigen::MatrixXd x(0, 3);
+  auto y = stan::math::row_map(
+      [](const Eigen::RowVectorXd& row) { return row; }, x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, row_map_empty_result) {
+  Eigen::MatrixXd x = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::row_map(
+      [&calls](const Eigen::RowVectorXd&) {
+        ++calls;
+        return Eigen::RowVectorXd(0);
+      },
+      x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, row_map_inconsistent_row_size) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  int call = 0;
+  EXPECT_THROW(stan::math::row_map(
+                   [&call](const Eigen::RowVectorXd&) {
+                     ++call;
+                     if (call == 1) {
+                       return Eigen::RowVectorXd::Ones(2);
+                     }
+                     return Eigen::RowVectorXd::Ones(3);
+                   },
+                   x),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, row_map_shared_arg_reused) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  auto offset = std::make_unique<double>(10.0);
+  auto y = stan::math::row_map(
+      [](const Eigen::RowVectorXd& row, const std::unique_ptr<double>& mu) {
+        return row.array() + *mu;
+      },
+      x, std::move(offset));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 12.0, 13.0, 14.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_map_sum_rows) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+
+  auto y = stan::math::col_map(
+      [](const Eigen::VectorXd& col) {
+        return Eigen::VectorXd::Constant(1, col.sum());
+      },
+      x);
+
+  Eigen::MatrixXd expected(1, 2);
+  expected << 4.0, 6.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_map_empty) {
+  Eigen::MatrixXd x(3, 0);
+  auto y
+      = stan::math::col_map([](const Eigen::VectorXd& col) { return col; }, x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, col_map_empty_result) {
+  Eigen::MatrixXd x = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::col_map(
+      [&calls](const Eigen::VectorXd&) {
+        ++calls;
+        return Eigen::VectorXd(0);
+      },
+      x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, col_map_inconsistent_col_size) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  int call = 0;
+  EXPECT_THROW(stan::math::col_map(
+                   [&call](const Eigen::VectorXd&) {
+                     ++call;
+                     if (call == 1) {
+                       return Eigen::VectorXd::Ones(2);
+                     }
+                     return Eigen::VectorXd::Ones(3);
+                   },
+                   x),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, col_map_shared_arg_reused) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  auto scale = std::make_unique<double>(2.0);
+  auto y = stan::math::col_map(
+      [](const Eigen::VectorXd& col, const std::unique_ptr<double>& s) {
+        return (*s) * col;
+      },
+      x, std::move(scale));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 2.0, 4.0, 6.0, 8.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_add) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+
+  auto y
+      = stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                                const Eigen::RowVectorXd& v) { return u + v; },
+                             a, b);
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_shared_arg_reused) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  auto offset = std::make_unique<double>(100.0);
+  auto y = stan::math::row_mapN(
+      [](const Eigen::RowVectorXd& u, const Eigen::RowVectorXd& v,
+         const std::unique_ptr<double>& mu) {
+        return *mu + u.array() + v.array();
+      },
+      std::forward_as_tuple(a, b), std::move(offset));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 111.0, 122.0, 133.0, 144.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_dim_mismatch) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 3);
+  b << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0;
+
+  EXPECT_THROW(
+      stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                              const Eigen::RowVectorXd& v) { return u + v; },
+                           a, b),
+      std::invalid_argument);
+}
+
+TEST(MathFunctor, row_mapN_empty) {
+  Eigen::MatrixXd a(0, 2);
+  Eigen::MatrixXd b(0, 2);
+  auto y
+      = stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                                const Eigen::RowVectorXd& v) { return u + v; },
+                             a, b);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, row_mapN_empty_result) {
+  Eigen::MatrixXd a = Eigen::MatrixXd::Ones(2, 2);
+  Eigen::MatrixXd b = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::row_mapN(
+      [&calls](const Eigen::RowVectorXd&, const Eigen::RowVectorXd&) {
+        ++calls;
+        return Eigen::RowVectorXd(0);
+      },
+      a, b);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, row_mapN_mixed_scalar_types) {
+  Eigen::MatrixXi a(2, 2);
+  a << 1, 2, 3, 4;
+  Eigen::MatrixXd b(2, 2);
+  b << 0.5, 1.5, 2.5, 3.5;
+  auto y = stan::math::row_mapN(
+      [](const auto& u, const auto& v) {
+        return u.template cast<double>() + v;
+      },
+      a, b);
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 1.5, 3.5, 5.5, 7.5;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_inconsistent_row_size) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  int call = 0;
+  EXPECT_THROW(
+      stan::math::row_mapN(
+          [&call](const Eigen::RowVectorXd&, const Eigen::RowVectorXd&) {
+            ++call;
+            if (call == 1) {
+              return Eigen::RowVectorXd::Ones(2);
+            }
+            return Eigen::RowVectorXd::Ones(3);
+          },
+          a, b),
+      std::invalid_argument);
+}
+
+TEST(MathFunctor, col_mapN_add) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v) { return u + v; },
+      a, b);
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_shared_arg_reused) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  auto offset = std::make_unique<double>(100.0);
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v,
+         const std::unique_ptr<double>& mu) {
+        return *mu + u.array() + v.array();
+      },
+      std::forward_as_tuple(a, b), std::move(offset));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 111.0, 122.0, 133.0, 144.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_dim_mismatch) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(3, 2);
+  b << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0;
+
+  EXPECT_THROW(
+      stan::math::col_mapN([](const Eigen::VectorXd& u,
+                              const Eigen::VectorXd& v) { return u + v; },
+                           a, b),
+      std::invalid_argument);
+}
+
+TEST(MathFunctor, col_mapN_empty) {
+  Eigen::MatrixXd a(2, 0);
+  Eigen::MatrixXd b(2, 0);
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v) { return u + v; },
+      a, b);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, col_mapN_empty_result) {
+  Eigen::MatrixXd a = Eigen::MatrixXd::Ones(2, 2);
+  Eigen::MatrixXd b = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::col_mapN(
+      [&calls](const Eigen::VectorXd&, const Eigen::VectorXd&) {
+        ++calls;
+        return Eigen::VectorXd(0);
+      },
+      a, b);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, col_mapN_mixed_scalar_types) {
+  Eigen::MatrixXi a(2, 2);
+  a << 1, 2, 3, 4;
+  Eigen::MatrixXd b(2, 2);
+  b << 0.5, 1.5, 2.5, 3.5;
+  auto y = stan::math::col_mapN(
+      [](const auto& u, const auto& v) {
+        return u.template cast<double>() + v;
+      },
+      a, b);
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 1.5, 3.5, 5.5, 7.5;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_inconsistent_col_size) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  int call = 0;
+  EXPECT_THROW(stan::math::col_mapN(
+                   [&call](const Eigen::VectorXd&, const Eigen::VectorXd&) {
+                     ++call;
+                     if (call == 1) {
+                       return Eigen::VectorXd::Ones(2);
+                     }
+                     return Eigen::VectorXd::Ones(3);
+                   },
+                   a, b),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, row_map_block_expression) {
+  Eigen::MatrixXd x(3, 3);
+  x << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0;
+  auto y = stan::math::row_map([](const auto& row) { return 2.0 * row; },
+                               x.block(0, 0, 2, 3));
+
+  Eigen::MatrixXd expected(2, 3);
+  expected << 2.0, 4.0, 6.0, 8.0, 10.0, 12.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_map_block_expression) {
+  Eigen::MatrixXd x(3, 3);
+  x << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0;
+  auto y = stan::math::col_map([](const auto& col) { return 2.0 * col; },
+                               x.block(0, 0, 3, 2));
+
+  Eigen::MatrixXd expected(3, 2);
+  expected << 2.0, 4.0, 8.0, 10.0, 14.0, 16.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_expression_inputs) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 0.0, 0.0, 1.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd c(2, 2);
+  c << 10.0, 20.0, 30.0, 40.0;
+
+  auto y
+      = stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                                const Eigen::RowVectorXd& v) { return u + v; },
+                             a * b, c);
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_expression_inputs) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 0.0, 0.0, 1.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd c(2, 2);
+  c << 10.0, 20.0, 30.0, 40.0;
+
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v) { return u + v; },
+      a * b, c);
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}