diff --git a/src/DCPSE/Dcpse.cuh b/src/DCPSE/Dcpse.cuh
index 327265041093a92c577c2b48b3c3e957aba47bda..fe3f2dc94f7834b1b1f7df24c9daba645aceb49c 100644
--- a/src/DCPSE/Dcpse.cuh
+++ b/src/DCPSE/Dcpse.cuh
@@ -49,6 +49,8 @@ public:
     // 2) The machinery for assembling and solving the linear system for coefficients starts...
     // 3) The user can then call an evaluate(point) method to get the evaluation of the differential operator
     //    on the given point.
+    ////c=HOverEpsilon. Note that the Eps value is computed by <h>/c (<h>=local average spacing for each particle and its support). This factor c is used in the Vandermonde.hpp.
+    double HOverEpsilon=0.9;
 private:
     const Point<dim, unsigned int> differentialSignature;
     const unsigned int differentialOrder;
@@ -501,7 +503,8 @@ private:
         const T condVTOL = 1e2;
 
         if (!isSharedSupport) {
-            SupportBuilder<vector_type> supportBuilder(particles, differentialSignature, rCut);
+            SupportBuilder<vector_type,vector_type>
+                supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
             unsigned int requiredSupportSize = monomialBasis.size();
             // need to resize supportKeys1D to yet unknown supportKeysTotalN
             // add() takes too long
@@ -524,12 +527,12 @@ private:
                 EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
                 // Vandermonde matrix computation
                 Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
-                        vandermonde(support, monomialBasis, particles);
+                        vandermonde(support, monomialBasis, particles, particles,HOverEpsilon);
+
                 vandermonde.getMatrix(V);
 
                 T condV = conditionNumber(V, condVTOL);
                 T eps = vandermonde.getEps();
-
                 if (condV > condVTOL) {
                     requiredSupportSize *= 2;
                     std::cout << "INFO: Increasing, requiredSupportSize = " << requiredSupportSize << std::endl; // debug
@@ -571,7 +574,8 @@ private:
         const T condVTOL = 1e2;
 
         if (!isSharedSupport) {
-            SupportBuilder<vector_type> supportBuilder(particles, differentialSignature, rCut);
+            SupportBuilder<vector_type,vector_type>
+                supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
             unsigned int requiredSupportSize = monomialBasis.size();
             // need to resize supportKeys1D to yet unknown supportKeysTotalN
             // add() takes too long
@@ -594,7 +598,8 @@ private:
                 EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
                 // Vandermonde matrix computation
                 Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
-                        vandermonde(support, monomialBasis, particles);
+                        vandermonde(support, monomialBasis, particles, particles, HOverEpsilon);
+
                 vandermonde.getMatrix(V);
 
                 T condV = conditionNumber(V, condVTOL);
@@ -660,7 +665,7 @@ std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution
                 size_t requiredSupportSize = monomialBasis.size() * supportSizeFactor;
                 // need to resize supportKeys1D to yet unknown supportKeysTotalN
                 // add() takes too long
-                SupportBuilder<vector_type> supportBuilder(particles, differentialSignature, rCut);
+                SupportBuilder<vector_type,vector_type> supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
                 kerOffsets.resize(supportRefs.size()+1);
 
                 while (it.isNext()) {
@@ -734,7 +739,7 @@ std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution
                 size_t requiredSupportSize = monomialBasis.size() * supportSizeFactor;
                 // need to resize supportKeys1D to yet unknown supportKeysTotalN
                 // add() takes too long
-                SupportBuilder<vector_type> supportBuilder(particles, differentialSignature, rCut);
+                SupportBuilder<vector_type,vector_type> supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
                 kerOffsets.resize(supportRefs.size()+1);
 
                 while (it.isNext()) {