diff --git a/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h b/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h
index e60b063107f..dc867e5b8de 100644
--- a/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h
+++ b/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h
@@ -124,12 +124,12 @@ namespace Heat_method_3 {
      * add `vd` to the source set, returning `false` if `vd` is already in the set.
      */
 
-     Matrix get_mass_matrix()
+     const Matrix& get_mass_matrix()
      {
        return mass_matrix;
      }
 
-     Matrix get_cotan_matrix()
+     const Matrix& get_cotan_matrix()
      {
        return cotan_matrix;
      }
@@ -158,7 +158,7 @@ namespace Heat_method_3 {
     /**
      *return current source set
     */
-    std::set<vertex_descriptor> get_sources()
+    const std::set<vertex_descriptor>& get_sources()
     {
       return sources;
     }
@@ -188,7 +188,7 @@ namespace Heat_method_3 {
     }
 
     /**
-     * get distance from the current source set to a vertex ` vd`.
+     * get distance from the current source set to a vertex `vd`.
      */
     double distance(vertex_descriptor vd)
     {
@@ -219,7 +219,7 @@ namespace Heat_method_3 {
     }
 
 
-    Matrix get_kronecker_delta()
+    const Matrix& get_kronecker_delta()
     {
       return kronecker;
     }
diff --git a/Heat_method_3/test/Heat_method_3/heat_method_surface_mesh_test.cpp b/Heat_method_3/test/Heat_method_3/heat_method_surface_mesh_test.cpp
index 3b3d8de5807..63ba30d10b9 100644
--- a/Heat_method_3/test/Heat_method_3/heat_method_surface_mesh_test.cpp
+++ b/Heat_method_3/test/Heat_method_3/heat_method_surface_mesh_test.cpp
@@ -43,8 +43,8 @@ int main()
   assert(!(hm.add_source(source)));
   assert(hm.remove_source(*(std::next(vertices(sm).first,3))));
   //cotan matrix tests
-  SparseMatrix M = hm.get_mass_matrix();
-  SparseMatrix c = hm.get_cotan_matrix();
+  const SparseMatrix& M = hm.get_mass_matrix();
+  const SparseMatrix& c = hm.get_cotan_matrix();
   double sum = 0;
   for(int k = 0; k<c.outerSize(); ++k)
   {
@@ -75,7 +75,7 @@ int main()
 
   //mass matrix tests
 
-  SparseMatrix K = hm.get_kronecker_delta();
+  const SparseMatrix& K = hm.get_kronecker_delta();
   assert(K.nonZeros()==1);