diff --git a/gtsam/base/Matrix.cpp b/gtsam/base/Matrix.cpp
index 9e56dfb6c..3cafdd0ba 100644
--- a/gtsam/base/Matrix.cpp
+++ b/gtsam/base/Matrix.cpp
@@ -184,21 +184,17 @@ void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, SubVec
 /* ************************************************************************* */
 //3 argument call
 void print(const Matrix& A, const string &s, ostream& stream) {
-  size_t m = A.rows(), n = A.cols();
-
-  // print out all elements
-  stream << s << "[\n";
-  for( size_t i = 0 ; i < m ; i++) {
-    for( size_t j = 0 ; j < n ; j++) {
-      double aij = A(i,j);
-      if(aij != 0.0)
-        stream << setw(12) << setprecision(9) << aij << ",\t";
-      else
-        stream << "         0.0,\t";
-    }
-    stream << endl;
-  }
-  stream << "];" << endl;
+  static const Eigen::IOFormat matlab(
+      Eigen::StreamPrecision, // precision
+      0, // flags
+      " ", // coeffSeparator
+      ";\n", // rowSeparator
+	  " \t",  // rowPrefix
+      "", // rowSuffix
+      "[\n", // matPrefix
+      "\n  ]" // matSuffix
+      );
+  cout << s << A.format(matlab) << endl;
 }
 
 /* ************************************************************************* */
diff --git a/gtsam/linear/JacobianFactor.cpp b/gtsam/linear/JacobianFactor.cpp
index 6e14960fd..d4df57298 100644
--- a/gtsam/linear/JacobianFactor.cpp
+++ b/gtsam/linear/JacobianFactor.cpp
@@ -347,13 +347,21 @@ JacobianFactor::JacobianFactor(const GaussianFactorGraph& graph,
 /* ************************************************************************* */
 void JacobianFactor::print(const string& s,
     const KeyFormatter& formatter) const {
+  static const Eigen::IOFormat matlab(
+      Eigen::StreamPrecision, // precision
+	  0, // flags
+      " ", // coeffSeparator
+      ";\n", // rowSeparator
+	  "\t",  // rowPrefix
+	  "", // rowSuffix
+      "[\n", // matPrefix
+      "\n  ]" // matSuffix
+      );
   if (!s.empty())
     cout << s << "\n";
   for (const_iterator key = begin(); key != end(); ++key) {
-    cout
-        << formatMatrixIndented(
-            (boost::format("  A[%1%] = ") % formatter(*key)).str(), getA(key))
-        << endl;
+    cout << boost::format("  A[%1%] = ") % formatter(*key);
+    cout << getA(key).format(matlab) << endl;
   }
   cout << formatMatrixIndented("  b = ", getb(), true) << "\n";
   if (model_)