diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt
index d90c68a317589123f8c0a0f1ad5af6d3b10a94c6..1cb974d741f108bbee0267d401e272b17f67c2cd 100644
--- a/unittest/CMakeLists.txt
+++ b/unittest/CMakeLists.txt
@@ -42,6 +42,7 @@ ADD_PYTHON_UNIT_TEST("py-matrix" "unittest/python/test_matrix.py" "unittest")
 ADD_PYTHON_UNIT_TEST("py-geometry" "unittest/python/test_geometry.py" "unittest")
 ADD_PYTHON_UNIT_TEST("py-complex" "unittest/python/test_complex.py" "unittest")
 ADD_PYTHON_UNIT_TEST("py-return-by-ref" "unittest/python/test_return_by_ref.py" "unittest")
+ADD_PYTHON_UNIT_TEST("py-eigen-ref" "unittest/python/test_eigen_ref.py" "unittest")
 
 ADD_PYTHON_UNIT_TEST("py-switch" "unittest/python/test_switch.py" "python/eigenpy")
 SET_TESTS_PROPERTIES("py-switch" PROPERTIES DEPENDS ${PYWRAP})
diff --git a/unittest/eigen_ref.cpp b/unittest/eigen_ref.cpp
index 32a972c195f29efc359ab75fe303b94358ea9282..1e58cfee4c7760169a605a79f3547b4a43369081 100644
--- a/unittest/eigen_ref.cpp
+++ b/unittest/eigen_ref.cpp
@@ -10,7 +10,7 @@ using namespace Eigen;
 using namespace eigenpy;
 
 template<typename MatType>
-void printMatrix(const Eigen::Ref<MatType> & mat)
+void printMatrix(const Eigen::Ref<const MatType> mat)
 {
   if(MatType::IsVectorAtCompileTime)
     std::cout << "isVector" << std::endl;
@@ -19,7 +19,7 @@ void printMatrix(const Eigen::Ref<MatType> & mat)
 }
 
 template<typename VecType>
-void printVector(const Eigen::Ref<VecType> & vec)
+void printVector(const Eigen::Ref<const VecType> & vec)
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VecType);
   printMatrix(vec);
@@ -28,9 +28,13 @@ void printVector(const Eigen::Ref<VecType> & vec)
 template<typename MatType>
 void setOnes(Eigen::Ref<MatType> mat)
 {
-  printMatrix(mat);
   mat.setOnes();
-  printMatrix(mat);
+}
+
+template<typename MatType>
+void fill(Eigen::Ref<MatType> mat, const typename MatType::Scalar & value)
+{
+  mat.fill(value);
 }
 
 BOOST_PYTHON_MODULE(eigen_ref)
@@ -48,4 +52,8 @@ BOOST_PYTHON_MODULE(eigen_ref)
   bp::def("setOnes", setOnes<Vector3d>);
   bp::def("setOnes", setOnes<VectorXd>);
   bp::def("setOnes", setOnes<MatrixXd>);
+  
+  bp::def("fillVec3", fill<Vector3d>);
+  bp::def("fillVec", fill<VectorXd>);
+  bp::def("fill", fill<MatrixXd>);
 }
diff --git a/unittest/python/test_eigen_ref.py b/unittest/python/test_eigen_ref.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe8cd29e386162e52b5292681d6611f228532524
--- /dev/null
+++ b/unittest/python/test_eigen_ref.py
@@ -0,0 +1,16 @@
+import numpy as np
+from eigen_ref import *
+
+def test(mat):
+
+  printMatrix(mat)
+  fill(mat,1.)
+  printMatrix(mat)
+  assert np.array_equal(mat,np.full(mat.shape,1.))
+
+rows = 10
+cols = 30
+
+mat = np.array(np.zeros((rows,cols)))
+
+test(mat)