diff --git a/Solver_interface/doc/Solver_interface/Solver_interface.txt b/Solver_interface/doc/Solver_interface/Solver_interface.txt
index 22fe8e4a5a7..ed50799101c 100644
--- a/Solver_interface/doc/Solver_interface/Solver_interface.txt
+++ b/Solver_interface/doc/Solver_interface/Solver_interface.txt
@@ -27,9 +27,7 @@ specific solver to be used must be given as template
 parameter.
 
 Each \cgal package using a sparse solver specifies which type of matrix
-and solver is required.
-
-Here are few examples: 
+and solver is required:
 
 \code{.cpp} 
 
@@ -46,6 +44,11 @@ typedef CGAL::Eigen_solver_traits< Eigen::SimplicialCholesky<EigenMatrix> > Dire
 
 \endcode
 
+Here is an example that shows how to fill the sparse matrix and call
+the solver:
+
+\cgalExample{Solver_interface/eigen_sparse_solvers.cpp}
+
 \section SectionSolverEigenSVD Eigen Singular Value Decomposition
 
 The class `Eigen_svd` provides an algorithm to solve in the least
diff --git a/Solver_interface/doc/Solver_interface/examples.txt b/Solver_interface/doc/Solver_interface/examples.txt
index a19fc9887e2..2d2b127bd9c 100644
--- a/Solver_interface/doc/Solver_interface/examples.txt
+++ b/Solver_interface/doc/Solver_interface/examples.txt
@@ -1,4 +1,5 @@
 /*!
 \example Solver_interface/eigen_singular_value_decomposition.cpp
 \example Solver_interface/eigen_variance_covariance_matrix.cpp
+\example Solver_interface/eigen_sparse_solvers.cpp
 */
diff --git a/Solver_interface/examples/Solver_interface/CMakeLists.txt b/Solver_interface/examples/Solver_interface/CMakeLists.txt
index b64849b0387..2d6de0e3734 100644
--- a/Solver_interface/examples/Solver_interface/CMakeLists.txt
+++ b/Solver_interface/examples/Solver_interface/CMakeLists.txt
@@ -26,6 +26,7 @@ if ( CGAL_FOUND )
 
   create_single_source_cgal_program( "eigen_singular_value_decomposition.cpp" )
   create_single_source_cgal_program( "eigen_variance_covariance_matrix.cpp" )
+  create_single_source_cgal_program( "eigen_sparse_solvers.cpp" )
 
 else()
   
diff --git a/Solver_interface/examples/Solver_interface/eigen_sparse_solvers.cpp b/Solver_interface/examples/Solver_interface/eigen_sparse_solvers.cpp
new file mode 100644
index 00000000000..4f26d5a284c
--- /dev/null
+++ b/Solver_interface/examples/Solver_interface/eigen_sparse_solvers.cpp
@@ -0,0 +1,64 @@
+#include <set>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Eigen_solver_traits.h>
+#include <CGAL/Eigen_matrix.h>
+
+
+typedef CGAL::Eigen_solver_traits<Eigen::ConjugateGradient<CGAL::Eigen_sparse_matrix<double>::EigenType> > Eigen_solver;
+typedef Eigen_solver::NT FT;
+typedef Eigen_solver::Matrix Eigen_matrix;
+typedef Eigen_solver::Vector Eigen_vector;
+
+
+int main(void)
+{
+  srand (time (NULL));
+  std::size_t degree = 3000;
+  std::size_t nb_nonzero_coef = 100;
+  
+  Eigen_vector B (degree);
+  Eigen_matrix A (degree);
+
+  std::set<std::pair<std::size_t, std::size_t> > done;
+  
+  Eigen_vector diag (degree);
+  for (std::size_t i = 0; i < nb_nonzero_coef; ++ i)
+    {
+      std::size_t x = rand () % degree;
+      std::size_t y = rand () % degree;
+      if (x == y)
+	continue;
+
+      FT value = rand () / static_cast<FT> (RAND_MAX);
+	
+      A.add_coef (x, y, value);
+      diag.set (x, diag.vector()[x] - value);
+      
+      B.set (x, 1.);
+    }
+
+  for (std::size_t i = 0; i < degree; ++ i)
+    A.add_coef (i, i, diag.vector()[i]);
+  
+  A.assemble_matrix();
+  
+  Eigen_vector X (degree);
+  FT d;
+
+  Eigen_solver solver;
+  if (!(solver.linear_solver (A, B, X, d)))
+    {
+      std::cerr << "Error: linear solver failed" << std::endl;
+      return -1;
+    }
+  
+  // Print extract of result
+  std::cout << "Vector X (non-zero coefficients) = [ ";
+  for (std::size_t i = 0; i < degree; ++ i)
+    if (X.vector()[i] != 0.)
+      std::cout << X.vector()[i] << " ";
+  std::cout << "]" << std::endl;
+  
+  return 0;
+}