diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/implicit_shape_smoothing_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/implicit_shape_smoothing_impl.h
index 78bce32d43c..920fc5794db 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/implicit_shape_smoothing_impl.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/implicit_shape_smoothing_impl.h
@@ -171,7 +171,6 @@ public:
   //void solve_system(Eigen_matrix& stiffness_matrix)
   void solve_system(Eigen::SparseMatrix<double>& L, const double time)
   {
-
     Eigen_matrix A(nb_vert_, nb_vert_);
 
     Eigen_vector bx(nb_vert_);
@@ -182,33 +181,33 @@ public:
     Eigen_vector Xy(nb_vert_);
     Eigen_vector Xz(nb_vert_);
 
-	std::cout << "compute coefficients...";
+	std::cerr << "compute coefficients...";
     compute_coeff_matrix(A, L, time);
-	std::cout << "done" << std::endl;
+	std::cerr << "done" << std::endl;
 
-	std::cout << "rhs...";
+	std::cerr << "rhs...";
 	compute_rhs(bx, by, bz);
-	std::cout << "done" << std::endl;
+	std::cerr << "done" << std::endl;
 
 
 	// Cholesky factorization
-	std::cout << "compute ...";
+	std::cerr << "compute ...";
 	Eigen_solver solver;
 	solver.compute(A);
-	std::cout << "done" << std::endl;
+	std::cerr << "done" << std::endl;
 	// solver.analyzePattern(A);
     // solver.factorize(A);
 
     // back-substitution
-	std::cout << "solve...";
+	std::cerr << "solve...";
 	Xx = solver.solve(bx);
     Xy = solver.solve(by);
     Xz = solver.solve(bz);
-	std::cout << "done" << std::endl;
+	std::cerr << "done" << std::endl;
 
     if(solver.info() != Eigen::Success)
     {
-      std::cerr<<"Not Solved!\n";
+      std::cerr << "Not Solved!" << std::endl;
       return;
     }
 
@@ -297,9 +296,13 @@ private:
       return C;
     }
 
-    Eigen_matrix stiff_matrix()
+    Eigen_matrix stiff_matrix() // TOFIX: pass reference to matrix
     {
-      Eigen_matrix mat(nb_vert_, nb_vert_);
+      Eigen_matrix mat(nb_vert_, nb_vert_); // TOFIX
+
+	  // TOFIX: very slow way to fill the matrx 
+	  // avoid matrix C!! and use triplets
+	  // read http://eigen.tuxfamily.org/dox/group__TutorialSparse.html#TutorialSparseFilling
 
       // cot values
       Eigen_matrix C = cot_entries();
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h
index db9b126f742..3de2bb0279c 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h
@@ -21,7 +21,10 @@ void smooth_shape(PolygonMesh& mesh, const double time, std::size_t nb_iter)
   internal::Shape_smoother<PolygonMesh, VertexPointMap> smoother(mesh, vpmap);
 
   Eigen::SparseMatrix<double> stiffness_matrix;
-  stiffness_matrix = smoother.calc_stiff_matrix();
+
+  std::cerr << "compute stiffness matrix...";
+  stiffness_matrix = smoother.calc_stiff_matrix(); // TOFIX: avoid copy
+  std::cerr << "done" << std::endl;
 
   for(std::size_t t=0; t<nb_iter; ++t)
   {