diff --git a/Surface_mesh_approximation/include/CGAL/VSA_approximation.h b/Surface_mesh_approximation/include/CGAL/VSA_approximation.h
index 204f31fd32f..d1e0ad5d3ef 100644
--- a/Surface_mesh_approximation/include/CGAL/VSA_approximation.h
+++ b/Surface_mesh_approximation/include/CGAL/VSA_approximation.h
@@ -518,7 +518,8 @@ public:
   }
 
   /*!
-   * @brief Merge two specified adjacent regions, the re-fitting is performed.
+   * @brief Merge two specified adjacent regions.
+   * The overall re-fitting is not performed and the proxy map is maintained.
    * @pre two proxies must be adjacent and px_enlarged < px_merged
    * @param px_enlarged the enlarged proxy
    * @param px_merged the merged proxy
@@ -617,26 +618,32 @@ public:
   }
 
   /*!
-   * @brief Merge adjacent regions, the re-fitting is performed.
-   * @param range_of_proxies range of proxies, must be adjacent
+   * @brief Split one proxy by default bisection, but N-section is also possible
+   * No re-fitting performed and the proxy map is maintained.
+   * @param px proxy index
+   * @param n split section
    * @return change of error
    */
-  // FT merge(range_of_proxies) {}
-  // Document what happens when the model has only one proxy,.
+  FT split(const std::size_t px, const std::size_t n = 2) {
+    if (px >= proxies.size())
+      return FT(0);
 
-  /*!
-   * @brief Split one proxy by default bisection, but N-section is also possible.
-   * @param p proxy
-   * @return change of error
-   */
-  // FT split(Proxy &p, int n = 2) {}
+    std::size_t count = 1;
+    FT err(0);
+    BOOST_FOREACH(face_descriptor f, faces(*m_pmesh)) {
+      if (count >= n)
+        break;
 
-  /*!
-   * @brief Split range of proxies by default bisection, but N-section is also possible.
-   * @param range_of_proxies range of proxies
-   * @return change of error
-   */
-  // FT split(range_of_proxies, int n = 2) {}
+      if (seg_pmap[f] == px && f != proxies[px].seed) {
+        err += (*fit_error)(f, proxies[px]);
+        seg_pmap[f] = proxies.size();
+        proxies.push_back(fit_new_proxy(f));
+        ++count;
+      }
+    }
+
+    return err;
+  }
 
   /*!
    * @brief Meshing, choose the default area weighted or the PCA plane fitting.
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
index 7d99da2618b..3f6109c1727 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
@@ -78,9 +78,20 @@ int main()
   // merge and teleport the proxies from local minimal
   std::cout << "teleport" << std::endl;
   l2_approx.teleport_proxies(3, false);
+  for (std::size_t i = 0; i < 10; ++i)
+    l2_approx.run_one_step();
   if (l2_approx.get_proxies_size() != 16)
     return EXIT_FAILURE;
 
+  // split proxy 0 into 2 proxies
+  // precondition: proxy 0 should have more than 2 facets
+  std::cout << "spliting" << std::endl;
+  l2_approx.split(0);
+  for (std::size_t i = 0; i < 10; ++i)
+    l2_approx.run_one_step();
+  if (l2_approx.get_proxies_size() != 17)
+    return EXIT_FAILURE;
+
   // extract the approximation polyhedron
   std::cout << "meshing" << std::endl;
   Polyhedron_3 out_mesh;