diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index dfd60694135c0c2cbb20f9ca84a17eb14ac4862d..f1b0cd940283596cdfb2f75254b4dabbbea0a04d 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -341,4 +341,20 @@ double linearSystemPETScBlockDouble::normInfRightHandSide() const
   return nor;
 }
 
+void linearSystemPETScBlockDouble::printMatlab(const char *filename) const
+{
+  _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
+  _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+  _try(VecAssemblyBegin(_b));
+  _try(VecAssemblyEnd(_b));
+
+  PetscViewer viewer;
+  PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer);
+  PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
+  printf("export mat to %s\n", filename);
+  MatView(_a, viewer);
+  PetscViewerDestroy(&viewer);
+  return;
+}
+
 #endif // HAVE_PETSC
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 040b062a50509e2b2989d94eae66dba8aeeb4faa..8dafd26108b663988a885be3055783e374c4ab1c 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -105,6 +105,7 @@ class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
   public:
   void _kspCreate();
   void print();
+  void printMatlab(const char *filename) const;
   virtual void addToMatrix(int row, int col, const fullMatrix<double> &val);
   virtual void addToRightHandSide(int row, const fullMatrix<double> &val);
   virtual void addToSolution(int row,  const fullMatrix<double> &val);