diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index f6cc81de812df86cdcbde1aa1a5d8669bf3ac428..c802a8e64b139ee38bc5584bc358c12ed3664796 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -220,14 +220,29 @@ class DCPSE_scheme: public MatrixAssembler
 
         // A and B must have the same rows
         if (row != row_b)
-            std::cerr << "Error " << __FILE__ << ":" << __LINE__ << "the term B and the Matrix A for Ax=B must contain the same number of rows\n";
+        {
+        	std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the term B and the Matrix A for Ax=B must contain the same number of rows\n";
+        	return;
+        }
+
+        if (row_b != p_map.size_local() * Sys_eqs::nvar)
+        {
+        	std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " your system is underdetermined you set " << row_b << " conditions " << " but i am expecting " << p_map.size_local() * Sys_eqs::nvar << std::endl;
+        	return;
+        }
 
         // Indicate all the non zero rows
         openfpm::vector<unsigned char> nz_rows;
         nz_rows.resize(row_b);
 
         for (size_t i = 0 ; i < trpl.size() ; i++)
+        {
+        	if (trpl.get(i).row() - s_pnt*Sys_eqs::nvar >= nz_rows.size())
+        	{
+        		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " It seems that you are setting colums that does not exist \n";
+        	}
             nz_rows.get(trpl.get(i).row() - s_pnt*Sys_eqs::nvar) = true;
+        }
 
         // Indicate all the non zero colums
         // This check can be done only on single processor
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index 34508295dbafa0fcea800cd79daa5874c1e84373..eedbace97392d03ec836dcf077d54394af022de1 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -176,12 +176,14 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
             Solver.impose(V_star[1], up_p,0,vy);
             Solver.impose(V_star[0], r_p, 0,vx);
             Solver.impose(V_star[1], r_p, 0,vy);
-            Solver.impose(V_star, dw_p,0,vx);
-            Solver.impose(V_star, l_p, 0,vy);
+            Solver.impose(V_star[0], dw_p,0,vx);
+            Solver.impose(V_star[1], dw_p,0,vy);
+            Solver.impose(V_star[0], l_p, 0,vx);
+            Solver.impose(V_star[1], l_p, 0,vy);
             Solver.solve(V_star);
             std::cout << "Stokes Solved" << std::endl;
             RHS=Div(V_star);
-            DCPSE_scheme<equations2,decltype(Particles)> SolverH(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+/*            DCPSE_scheme<equations2,decltype(Particles)> SolverH(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
             auto Helmholtz = Lap(H);
             auto D_y=Dy(H);
             auto D_x=Dx(H);
@@ -194,7 +196,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
             std::cout << "Helmholtz Solved" << std::endl;
             V=V_star-Grad(H);
             P=P+Lap(H);
-            std::cout << "V,P Corrected" << std::endl;
+            std::cout << "V,P Corrected" << std::endl;*/
             Particles.write_frame("Stokes",i);
 
         }