diff --git a/src/FiniteDifference/Derivative.hpp b/src/FiniteDifference/Derivative.hpp
index d0df1be44a843891fa688ff25c373c370d86026d..240d6507ba411d08fcb2199cb9b794e45620bebc 100644
--- a/src/FiniteDifference/Derivative.hpp
+++ b/src/FiniteDifference/Derivative.hpp
@@ -31,7 +31,7 @@ class D
 	 * \tparam ord
 	 *
 	 */
-	inline static std::unordered_map<long int,typename Sys_eqs::stype> value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+	inline static void value(const grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
 	{
 		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
 	}
diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index ea0cc0e7ad2b0d4b27b3476a9439017e8756ed49..ea5154bd61517dd6661e53133616e65a35fb3bad 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -46,7 +46,7 @@ public:
 			// convert into global coordinate the position
 			auto keyg = it.getGKey(key);
 
-			// Convert local key into global key
+			// Calculate the non-zero colums
 			T::value(keyg,gs,cols,1.0);
 
 			// create the triplet
diff --git a/src/FiniteDifference/eq.hpp b/src/FiniteDifference/eq.hpp
index ea83473e6bcfa306854afc5cdda926e3d7c19e74..321d0c39f3958ae6a535d054e16e03364c867777 100644
--- a/src/FiniteDifference/eq.hpp
+++ b/src/FiniteDifference/eq.hpp
@@ -103,49 +103,6 @@ class minus
 	}
 };
 
-/*! \brief It model an expression expr1 * expr2
- *
- * \warning expr1 MUST be a constant expression
- *
- * \tparam expr1
- * \tparam expr2
- *
- */
-template<typename expr1,typename expr2, typename Sys_eqs>
-class mul
-{
-	/*! \brief Create the row of the Matrix
-	 *
-	 * \tparam ord
-	 *
-	 */
-	template<unsigned int ord=EQS_FIELD> static void value(const grid_key_dx<Sys_eqs::dims> & pos)
-	{
-		if (EQS_FIELD)
-			value_f(pos);
-		else
-			value_s(pos);
-	}
-
-	/*! \brief fill the row
-	 *
-	 *
-	 */
-	static openfpm::vector<triplet<typename Sys_eqs::stype>> value_s(grid_key_dx<Sys_eqs::dims> & it)
-	{
-		return expr1::const_value(it) * expr2::value_s(it);
-	}
-
-	/*! \brief fill the row
-	 *
-	 *
-	 */
-	static void value_f(grid_key_dx<Sys_eqs::dims> & it)
-	{
-		return expr1::const_value(it) * expr2::value_s(it);
-	}
-};
-
 // spatial position + value
 
 template<unsigned int dim,typename T>
diff --git a/src/FiniteDifference/eq_unit_test.hpp b/src/FiniteDifference/eq_unit_test.hpp
index e9f97edaefb4e5a824ed70e625f17ffe9d2d351c..4a6fb1ff8aa9dd490a285b1c9c8487c3a2d5521f 100644
--- a/src/FiniteDifference/eq_unit_test.hpp
+++ b/src/FiniteDifference/eq_unit_test.hpp
@@ -11,6 +11,7 @@
 #include "Laplacian.hpp"
 #include "FiniteDifference/eq.hpp"
 #include "FiniteDifference/sum.hpp"
+#include "FiniteDifference/mul.hpp"
 #include "Grid/grid_dist_id.hpp"
 #include "data_type/scalar.hpp"
 #include "Decomposition/CartDecomposition.hpp"
@@ -41,12 +42,12 @@ struct lid_nn
 const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
 
 // Constant Field
-
 struct eta
 {
+	float val()	{return 1.0;}
 };
 
-// Model the equation
+// Model the equations
 
 constexpr unsigned int v[] = {0,1};
 constexpr unsigned int P = 2;
@@ -73,6 +74,8 @@ typedef D<x,v_x,lid_nn> dx_vx;
 typedef D<y,v_y,lid_nn> dy_vy;
 typedef sum<dx_vx,dy_vy> incompressibility;
 
+BOOST_AUTO_TEST_SUITE( eq_test_suite )
+
 // Lid driven cavity, uncompressible fluid
 
 BOOST_AUTO_TEST_CASE( lid_driven_cavity )
@@ -107,4 +110,6 @@ BOOST_AUTO_TEST_CASE( lid_driven_cavity )
 	fd.impose(vy, g_dist.getGridInfo(), g_dist.getSubDomainIterator(bulk_start,bulk_end));
 }
 
+BOOST_AUTO_TEST_SUITE_END()
+
 #endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_ */
diff --git a/src/FiniteDifference/sum.hpp b/src/FiniteDifference/sum.hpp
index f2eca613e5c21e81415f655e3d82ab338d608cb0..0f4eb14815de3a949df10382ae9b2f3ae2a8366d 100644
--- a/src/FiniteDifference/sum.hpp
+++ b/src/FiniteDifference/sum.hpp
@@ -9,6 +9,46 @@
 #define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_SUM_HPP_
 
 #include <boost/mpl/vector.hpp>
+#include "config.h"
+#include <unordered_map>
+#include "util/for_each_ref.hpp"
+
+template<typename v_expr>
+struct sum_functor_value
+{
+	//! Number of elements in the vector v_expr
+	typedef boost::mpl::size<v_expr> size;
+
+	//! Last element of sum
+	typedef typename boost::mpl::at<v_expr,boost::mpl::int_<size::value-1> >::type last;
+
+	//! sum functor
+	std::unordered_map<long int,typename last::stype> & cols;
+
+	const grid_sm<last::dims,void> & gs;
+
+	//! position
+	grid_key_dx<last::dims> & key;
+
+	//! coefficent
+	typename last::stype coeff;
+
+	/*! \brief constructor
+	 *
+	 */
+	sum_functor_value(grid_key_dx<last::dims> & key, const grid_sm<last::dims,void> & gs, typename last::stype coeff)
+	:cols(cols),gs(gs),key(key),coeff(coeff)
+	{};
+
+
+
+	//! It call this function for every expression in the sum
+	template<typename T>
+	void operator()(T& t) const
+	{
+		boost::mpl::at<v_expr, boost::mpl::int_<T::value> >::type::value(key,gs,cols,coeff);
+	}
+};
 
 /*! \brief It model an expression expr1 + ... exprn
  *
@@ -32,9 +72,13 @@ struct sum
 	 * \tparam ord
 	 *
 	 */
-	inline static std::unordered_map<long int,typename Sys_eqs::stype> value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+	inline static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
 	{
-		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
+		// Sum functor
+		sum_functor_value<v_expr> sm(pos,gs,coeff);
+
+		// for each element in the expression calculate the non-zero Matrix elements
+		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,v_sz::type::value - 1> >(sm);
 	}
 
 	/*! \brief Calculate the position where the derivative is calculated
diff --git a/src/main.cpp b/src/main.cpp
index 26c8feef9e8b06bd082707ceecfb83bb0ddc7ea8..e8e3e392f5b8d52eca1399dcb56754ae53a3f1b4 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,3 +6,4 @@
 #include "FiniteDifference/FDScheme_unit_tests.hpp"
 #include "FiniteDifference/util/common_test.hpp"
 #include "FiniteDifference/eq_unit_test.hpp"
+#include "util/util_num_unit_tests.hpp"