diff --git a/include/cantera/numerics/DenseMatrix.h b/include/cantera/numerics/DenseMatrix.h
index 7f4b2d6c4..2ccea46a8 100644
--- a/include/cantera/numerics/DenseMatrix.h
+++ b/include/cantera/numerics/DenseMatrix.h
@@ -112,8 +112,9 @@ public:
 
     //! Multiply A*B and write result to \c prod.
     /*!
-     * @param[in]  b     DenseMatrix B of size NxN
-     * @param[out] prod  DenseMatrix prod size NxN
+     * Take this matrix to be of size NxM.
+     * @param[in]  b     DenseMatrix B of size MxP
+     * @param[out] prod  DenseMatrix prod size NxP
      */
     virtual void mult(const DenseMatrix& b, DenseMatrix& prod) const;
 
diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp
index a716e349b..53b84c34c 100644
--- a/src/numerics/DenseMatrix.cpp
+++ b/src/numerics/DenseMatrix.cpp
@@ -93,13 +93,18 @@ void DenseMatrix::mult(const double* b, double* prod) const
 
 void DenseMatrix::mult(const DenseMatrix& B, DenseMatrix& prod) const
 {
-    if (m_ncols != B.nColumns() || m_nrows != B.nRows() || m_ncols != m_nrows || m_ncols != prod.nColumns() || m_nrows != prod.nColumns()) {
-        throw CanteraError("mult(const DenseMatrix &B, DenseMatrix &prod)",
-                           "Cannot multiply matrices that are not square and/or not the same size.");
+    if (nColumns() != B.nRows()) {
+        throw CanteraError("DenseMatrix::mult",
+            "Inner matrix dimensions do not agree: {} != {}", nColumns(), B.nRows());
+    }
+    if (nRows() != prod.nRows() || B.nColumns() != prod.nColumns()) {
+        throw CanteraError("DenseMatrix::mult",
+            "Output matrix has wrong dimensions: {}x{} != {}x{}",
+            prod.nRows(), prod.nColumns(), nRows(), B.nColumns());
     }
     const doublereal* const* bcols = B.const_colPts();
     doublereal* const* prodcols = prod.colPts();
-    for (size_t col=0; col < m_ncols; ++col) {
+    for (size_t col=0; col < B.nColumns(); ++col) {
         // Loop over ncols multiplying A*column of B and storing in
         // corresponding prod column
         mult(bcols[col], prodcols[col]);
diff --git a/test/general/test_matrices.cpp b/test/general/test_matrices.cpp
index b7a140a7f..b4caa747c 100644
--- a/test/general/test_matrices.cpp
+++ b/test/general/test_matrices.cpp
@@ -129,6 +129,16 @@ public:
         A1(3,3) = -9;
     }
 
+    double special_sum(DenseMatrix M) {
+        double sum = 0;
+        for (size_t i = 0; i < M.nRows(); i++) {
+            for (size_t j = 0; j < M.nColumns(); j++) {
+                sum += M(i,j) * (i+2*j+1);
+            }
+        }
+        return sum;
+    }
+
     DenseMatrix A1, A2, A3;
     vector_fp x4, x3;
     vector_fp b1, b2, b3;
@@ -150,3 +160,22 @@ TEST_F(DenseMatrixTest, matrix_times_vector)
         EXPECT_DOUBLE_EQ(b3[i], c[i]);
     }
 }
+
+TEST_F(DenseMatrixTest, matrix_times_matrix)
+{
+    DenseMatrix c(4, 3);
+    A1.mult(A2, c);
+    EXPECT_DOUBLE_EQ(3033, special_sum(c));
+
+    c.resize(3, 4);
+    A3.mult(A1, c);
+    EXPECT_DOUBLE_EQ(3386, special_sum(c));
+
+    c.resize(3, 3);
+    A3.mult(A2, c);
+    EXPECT_DOUBLE_EQ(2989, special_sum(c));
+
+    c.resize(4, 4);
+    A2.mult(A3, c);
+    EXPECT_DOUBLE_EQ(4014, special_sum(c));
+}