Index: tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.Grid.tol
===================================================================
--- tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.Grid.tol	(revision 38052b371de6c0ad00cfa69bfbb6db3ba1b33364)
+++ tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.Grid.tol	(revision 3f198ccdfd9b11b0fcc7996a8a441752bf89b410)
@@ -14,12 +14,17 @@
 //////////////////////////////////////////////////////////////////////////////
 {
+  @IndexedGrid.Border _.grid;
 //////////////////////////////////////////////////////////////////////////////
 //Configuration parameters
 //////////////////////////////////////////////////////////////////////////////
-  @IndexedGrid.Border _.grid;
+  Real _.radius;
+  Real _.radius.margin;
+  Real _.grid.size;
 
 //////////////////////////////////////////////////////////////////////////////
 //Auxiliar members
 //////////////////////////////////////////////////////////////////////////////
+  VMatrix _.Y.h = Rand(0,0,0,0);
+  VMatrix _.Y.neighbourhood = Rand(0,0,0,0);
 
   ////////////////////////////////////////////////////////////////////////////
@@ -45,4 +50,11 @@
 
   //////////////////////////////////////////////////////////////////////////////
+  VMatrix get.radius(VMatrix Y.idx)
+  //////////////////////////////////////////////////////////////////////////////
+  {
+    Constant(VRows(Y.idx),1,_.radius)
+  };
+
+  //////////////////////////////////////////////////////////////////////////////
   VMatrix closest.reference.index(VMatrix X)
   //////////////////////////////////////////////////////////////////////////////
@@ -68,4 +80,79 @@
     VMatrix _.Yj := Mat2VMat(SetMat([[EvalSet(X.class,Card)]]),1);
     Real _.H := _.grid::_.H
+  };
+
+  //////////////////////////////////////////////////////////////////////////////
+  Real build.stress.points(Real void)
+  //////////////////////////////////////////////////////////////////////////////
+  {
+    Text _MID = "[@DesignC2.Spheric::build.stress.points] ";
+    WriteLn("TRACE "+_MID+"1");
+    Matrix Y.neighbour.triplet = Group("ConcatRows", For(1,_.H,Matrix(Real h)
+    {
+      Set neighbours = Unique(VRow2Set(SubRow(_.Y.neighbourhood,[[h]])));
+      Set neighbours.up = Select(neighbours,Real(Real s) { s>h });
+      Matrix s = SetCol(Sort(neighbours.up,Compare));
+      Rand(Rows(s),1,h,h)|s|Rand(Rows(s),1,1,1)
+    }));
+
+    WriteLn("TRACE "+_MID+"2");
+    VMatrix Y.neighbour.1.map = Triplet(Y.neighbour.triplet,_.H,_.H);
+    VMatrix Y.neighbour.2.map = Y.neighbour.1.map*Y.neighbour.1.map;
+    VMatrix Y.neighbour.12.map = Max(Y.neighbour.1.map,Y.neighbour.2.map);
+    Matrix Y.neighbour.12.triplet = VMat2Triplet(Y.neighbour.12.map);
+    VMatrix Z.h.a = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[1]]));
+    VMatrix Z.h.b = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[2]]));
+    VMatrix Y.neighbour.a = SubRow(_.Y,VCol2Set(Z.h.a));
+    VMatrix Y.neighbour.b = SubRow(_.Y,VCol2Set(Z.h.b));
+    VMatrix Z.all = (Y.neighbour.a + Y.neighbour.b)*0.5;
+
+    WriteLn("TRACE "+_MID+"3");
+    VMatrix Z.interior.a = 
+    {
+      VMatrix Z.norm2 = RPow(Y.neighbour.a-Z.all,2)*Rand(_.n,1,1,1);
+      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
+      LE(Z.tao,Rand(VRows(Z.tao),1,1,1)) 
+    };
+    VMatrix Z.interior.b = 
+    {
+      VMatrix Z.norm2 = RPow(Y.neighbour.b-Z.all,2)*Rand(_.n,1,1,1);
+      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
+      LE(Z.tao,Rand(VRows(Z.tao),1,1,1)) 
+    };
+    VMatrix match.interior = Z.interior.a $* Z.interior.b;
+
+    Set sel.interior = MatQuery::SelectMatch(match.interior);
+    VMatrix Z.interior = SubRow(Z.all,sel.interior);
+    VMatrix Z.h.interior = SubRow(Z.h.a | Z.h.b,sel.interior);
+    Set criterium = Range(1,_.n,1);
+    Set sel.unique = Col2Set(Sort(VMat2Mat(Z.interior),criterium,True));
+    WriteLn("TRACE "+_MID+"4");
+    VMatrix _.Z := SubRow(Z.interior, sel.unique);
+    VMatrix _.Z.h := SubRow(Z.h.interior, sel.unique);
+    Real _.Tz := VRows(_.Z);
+    VMatrix _.Z.Y.neighbour := 
+    {
+      Matrix Z.t = DifEq(1/(1-B),Rand(_.Tz,1,1,1));
+      Matrix Z.th = (Z.t << Z.t) | 
+       VMat2Mat((SubCol(_.Z.h,[[1]])<<SubCol(_.Z.h,[[2]])));
+      Matrix Z.Y.idx = PivotByRows(Z.th,Sort(Z.th,[[2]]));
+      VMatrix Z.Y.dif = 
+       SubRow(_.Z,Col2Set(SubCol(Z.Y.idx,[[1]]))) -
+       SubRow(_.Y,Col2Set(SubCol(Z.Y.idx,[[2]])));
+        VMatrix Z.Y.norm2 = RPow(Z.Y.dif,2)*Rand(_.n,1,1,1);
+        VMatrix Z.Y.tao = RPow(Z.Y.norm2,0.5)/_.radius;
+        VMatrix Z.Y.match = LE(Z.Y.tao,Rand(VRows(Z.Y.tao),1,1,1));
+      SubRow(Mat2VMat(Z.Y.idx),MatQuery::SelectMatch(Z.Y.match))
+    };
+    VMatrix aux = MatQuery::SortAndCount(SubCol(_.Z.Y.neighbour,[[2]]));
+    VMatrix _.Yt := Triplet(VMat2Mat(SubCol(aux,[[1]])) |  
+                            Rand(VRows(aux),1,1,1) | 
+                            VMat2Mat(SubCol(aux,[[2]])),_.H,1);
+    Real _.Mz := Case(
+      prior.Z.order< 0, 0,
+      prior.Z.order==0, _.Tz,
+      prior.Z.order==1, _.Tz*(1+_.n),
+      prior.Z.order==2, _.Tz*(1+_.n*3/2+(1/2)*_.n^2));
+    _.Tz
   };
 
Index: tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.MeshFree.tol
===================================================================
--- tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.MeshFree.tol	(revision 38052b371de6c0ad00cfa69bfbb6db3ba1b33364)
+++ tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.MeshFree.tol	(revision 3f198ccdfd9b11b0fcc7996a8a441752bf89b410)
@@ -2,5 +2,5 @@
 //////////////////////////////////////////////////////////////////////////////
 // FILE   : DesignC2.Monomial.Spheric.MeshFree.tol
-// PURPOSE: Defines Class @DesignC2.Monomial.Spheric.MeshFree.tol
+// PURPOSE: Defines Class @DesignC2.Monomial.Spheric.MeshFree
 //////////////////////////////////////////////////////////////////////////////
 
@@ -9,5 +9,5 @@
 Text _.autodoc.member.DesignC2.Monomial.Spheric.MeshFree.tol = 
 "";
-Class @DesignC2.Monomial.Spheric.MeshFree.tol : 
+Class @DesignC2.Monomial.Spheric.MeshFree : 
   @DesignC2.Monomial, 
   @DesignC2.Spheric
@@ -17,5 +17,5 @@
 //Configuration parameters
 //////////////////////////////////////////////////////////////////////////////
-  @IndexedGrid.Border _.grid;
+  MatQuery::@ANN.KDTree _.ann;
 
 //////////////////////////////////////////////////////////////////////////////
@@ -24,5 +24,5 @@
 
   ////////////////////////////////////////////////////////////////////////////
-  Static @DesignC2.Monomial.Spheric.MeshFree.tol New(
+  Static @DesignC2.Monomial.Spheric.MeshFree New(
     VMatrix X, 
     Real radius.margin,
@@ -31,5 +31,5 @@
   ////////////////////////////////////////////////////////////////////////////
   {
-    @DesignC2.Monomial.Spheric.MeshFree.tol new =
+    @DesignC2.Monomial.Spheric.MeshFree new =
     [[ 
       VMatrix _.X = X;
@@ -78,5 +78,5 @@
       Real _.Tb := 0
     },{
-    Text _MID = "[@DesignC2.Monomial.Spheric.MeshFree.tol::build.border.points] ";
+    Text _MID = "[@DesignC2.Monomial.Spheric.MeshFree::build.border.points] ";
     VMatrix   B.t  = Mat2VMat(SubCol(_.grid::_.border.assign,[[1]]));
     VMatrix _.B.h := Mat2VMat(SubCol(_.grid::_.border.assign,[[2]]));
@@ -101,5 +101,5 @@
   //////////////////////////////////////////////////////////////////////////////
   {
-    Text _MID = "[@DesignC2.Monomial.Spheric.MeshFree.tol::build.model.design] ";
+    Text _MID = "[@DesignC2.Monomial.Spheric.MeshFree::build.model.design] ";
     FunRn2R::@BaseC2 base.max = 
       FunRn2R::@BaseC2.Monomial::New(_.n,max.deg);
Index: tol_pkg/MulVarFunAprox/DesignC2.Spheric.tol
===================================================================
--- tol_pkg/MulVarFunAprox/DesignC2.Spheric.tol	(revision 38052b371de6c0ad00cfa69bfbb6db3ba1b33364)
+++ tol_pkg/MulVarFunAprox/DesignC2.Spheric.tol	(revision 3f198ccdfd9b11b0fcc7996a8a441752bf89b410)
@@ -12,94 +12,11 @@
 //////////////////////////////////////////////////////////////////////////////
 {
-//////////////////////////////////////////////////////////////////////////////
-//Configuration parameters
-//////////////////////////////////////////////////////////////////////////////
-  Real _.radius;
-  Real _.radius.margin;
-  Real _.grid.size;
-
-//////////////////////////////////////////////////////////////////////////////
-//Auxiliar members
-//////////////////////////////////////////////////////////////////////////////
-  VMatrix _.Y.h = Rand(0,0,0,0);
-  VMatrix _.Y.neighbourhood = Rand(0,0,0,0);
-
   //////////////////////////////////////////////////////////////////////////////
-  Real build.stress.points(Real void)
+  //Virtual methods
   //////////////////////////////////////////////////////////////////////////////
-  {
-    Text _MID = "[@DesignC2.Spheric::build.stress.points] ";
-    WriteLn("TRACE "+_MID+"1");
-    Matrix Y.neighbour.triplet = Group("ConcatRows", For(1,_.H,Matrix(Real h)
-    {
-      Set neighbours = Unique(VRow2Set(SubRow(_.Y.neighbourhood,[[h]])));
-      Set neighbours.up = Select(neighbours,Real(Real s) { s>h });
-      Matrix s = SetCol(Sort(neighbours.up,Compare));
-      Rand(Rows(s),1,h,h)|s|Rand(Rows(s),1,1,1)
-    }));
-
-    WriteLn("TRACE "+_MID+"2");
-    VMatrix Y.neighbour.1.map = Triplet(Y.neighbour.triplet,_.H,_.H);
-    VMatrix Y.neighbour.2.map = Y.neighbour.1.map*Y.neighbour.1.map;
-    VMatrix Y.neighbour.12.map = Max(Y.neighbour.1.map,Y.neighbour.2.map);
-    Matrix Y.neighbour.12.triplet = VMat2Triplet(Y.neighbour.12.map);
-    VMatrix Z.h.a = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[1]]));
-    VMatrix Z.h.b = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[2]]));
-    VMatrix Y.neighbour.a = SubRow(_.Y,VCol2Set(Z.h.a));
-    VMatrix Y.neighbour.b = SubRow(_.Y,VCol2Set(Z.h.b));
-    VMatrix Z.all = (Y.neighbour.a + Y.neighbour.b)*0.5;
-
-    WriteLn("TRACE "+_MID+"3");
-    VMatrix Z.interior.a = 
-    {
-      VMatrix Z.norm2 = RPow(Y.neighbour.a-Z.all,2)*Rand(_.n,1,1,1);
-      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
-      LE(Z.tao,Rand(VRows(Z.tao),1,1,1)) 
-    };
-    VMatrix Z.interior.b = 
-    {
-      VMatrix Z.norm2 = RPow(Y.neighbour.b-Z.all,2)*Rand(_.n,1,1,1);
-      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
-      LE(Z.tao,Rand(VRows(Z.tao),1,1,1)) 
-    };
-    VMatrix match.interior = Z.interior.a $* Z.interior.b;
-
-    Set sel.interior = MatQuery::SelectMatch(match.interior);
-    VMatrix Z.interior = SubRow(Z.all,sel.interior);
-    VMatrix Z.h.interior = SubRow(Z.h.a | Z.h.b,sel.interior);
-    Set criterium = Range(1,_.n,1);
-    Set sel.unique = Col2Set(Sort(VMat2Mat(Z.interior),criterium,True));
-    WriteLn("TRACE "+_MID+"4");
-    VMatrix _.Z := SubRow(Z.interior, sel.unique);
-    VMatrix _.Z.h := SubRow(Z.h.interior, sel.unique);
-    Real _.Tz := VRows(_.Z);
-    VMatrix _.Z.Y.neighbour := 
-    {
-      Matrix Z.t = DifEq(1/(1-B),Rand(_.Tz,1,1,1));
-      Matrix Z.th = (Z.t << Z.t) | 
-       VMat2Mat((SubCol(_.Z.h,[[1]])<<SubCol(_.Z.h,[[2]])));
-      Matrix Z.Y.idx = PivotByRows(Z.th,Sort(Z.th,[[2]]));
-      VMatrix Z.Y.dif = 
-       SubRow(_.Z,Col2Set(SubCol(Z.Y.idx,[[1]]))) -
-       SubRow(_.Y,Col2Set(SubCol(Z.Y.idx,[[2]])));
-        VMatrix Z.Y.norm2 = RPow(Z.Y.dif,2)*Rand(_.n,1,1,1);
-        VMatrix Z.Y.tao = RPow(Z.Y.norm2,0.5)/_.radius;
-        VMatrix Z.Y.match = LE(Z.Y.tao,Rand(VRows(Z.Y.tao),1,1,1));
-      SubRow(Mat2VMat(Z.Y.idx),MatQuery::SelectMatch(Z.Y.match))
-    };
-    VMatrix aux = MatQuery::SortAndCount(SubCol(_.Z.Y.neighbour,[[2]]));
-    VMatrix _.Yt := Triplet(VMat2Mat(SubCol(aux,[[1]])) |  
-                            Rand(VRows(aux),1,1,1) | 
-                            VMat2Mat(SubCol(aux,[[2]])),_.H,1);
-    Real _.Mz := Case(
-      prior.Z.order< 0, 0,
-      prior.Z.order==0, _.Tz,
-      prior.Z.order==1, _.Tz*(1+_.n),
-      prior.Z.order==2, _.Tz*(1+_.n*3/2+(1/2)*_.n^2));
-    _.Tz
-  };
+  VMatrix get.radius(VMatrix Y.idx);
 
   ///////////////////////////////////////////////////////////////////////////
-  Set get.X.weight(Set X.dif)
+  Set get.X.weight(Set X.dif, Set X.h)
   ///////////////////////////////////////////////////////////////////////////
   {
@@ -108,5 +25,5 @@
     {
       VMatrix X.norm2 = RPow(X.dif[i],2)*Rand(_.n,1,1,1);
-      RPow(X.norm2,0.5)/_.radius
+      RPow(X.norm2,0.5)$/get.radius(X.h[i])
     });
     Set X.tao2 = For(1,_.X.maxNeighbour,VMatrix(Real i)
@@ -130,5 +47,5 @@
   };
   ///////////////////////////////////////////////////////////////////////////
-  Set get.X.weight.partial.derivative.1(Set X.dif, Real s)
+  Set get.X.weight.partial.derivative.1(Set X.dif, Set X.h, Real s)
   ///////////////////////////////////////////////////////////////////////////
   {
@@ -137,5 +54,5 @@
     {
       VMatrix X.norm2 = RPow(X.dif[i],2)*Rand(_.n,1,1,1);
-      RPow(X.norm2,0.5)/_.radius
+      RPow(X.norm2,0.5)$/get.radius(X.h[i])
     });
     Set X.tao2 = For(1,_.X.maxNeighbour,VMatrix(Real i)
@@ -153,5 +70,7 @@
     Set X.lambda.der = For(1,_.X.maxNeighbour,VMatrix(Real i)
     {
-      RPow(X.match[i] $* (Rand(m,1,1,1)-X.tao2[i]),3)$*SubCol(X.dif[i],[[s]])*8/_.radius
+      RPow(X.match[i] $* 
+           (Rand(m,1,1,1)-X.tao2[i]),3) $*
+           (SubCol(X.dif[i],[[s]])*8) $/ get.radius(X.h[i])
     });
     VMatrix X.lambda.sum = SetSum(X.lambda);
Index: tol_pkg/MulVarFunAprox/DesignC2.tol
===================================================================
--- tol_pkg/MulVarFunAprox/DesignC2.tol	(revision 38052b371de6c0ad00cfa69bfbb6db3ba1b33364)
+++ tol_pkg/MulVarFunAprox/DesignC2.tol	(revision 3f198ccdfd9b11b0fcc7996a8a441752bf89b410)
@@ -75,6 +75,6 @@
   VMatrix closest.reference.index(VMatrix X);
   VMatrix match.evaluable(VMatrix X);
-  Set get.X.weight(Set X.dif);
-  Set get.X.weight.partial.derivative.1(Set X.dif, Real s);
+  Set get.X.weight(Set X.dif, Set X.h);
+  Set get.X.weight.partial.derivative.1(Set X.dif, Set X.h, Real s);
 
   //////////////////////////////////////////////////////////////////////////////
@@ -118,11 +118,12 @@
     VMatrix X.j=DifEq(1/(1-B),Rand(m,1,1,1));
     WriteLn("TRACE "+_MID+"3");
+    Set X.h_ = For(1,_.X.maxNeighbour,VMatrix(Real i) { SubCol(X.h,[[i]]) });
     Set X.dif = For(1,_.X.maxNeighbour,VMatrix(Real i)
     {
-      Set X.h.i = VCol2Set(SubCol(X.h,[[i]]));
+      Set X.h.i = VCol2Set(X.h_[i]);
       X - SubRow(_.Y,X.h.i)
     });
     WriteLn("TRACE "+_MID+"4");
-    Set X.weight = get.X.weight(X.dif);
+    Set X.weight = get.X.weight(X.dif, X.h_);
     WriteLn("TRACE "+_MID+"5");
     VMatrix YnCum = DifEq(1/(1-B),_.Yn);
@@ -169,12 +170,13 @@
     VMatrix X.j=DifEq(1/(1-B),Rand(m,1,1,1));
     WriteLn("TRACE "+_MID+"4");
+    Set X.h_ = For(1,_.X.maxNeighbour,VMatrix(Real i) { SubCol(X.h,[[i]]) });
     Set X.dif = For(1,_.X.maxNeighbour,VMatrix(Real i)
     {
-      Set X.h.i = VCol2Set(SubCol(X.h,[[i]]));
+      Set X.h.i = VCol2Set(X.h_[i]);
       X - SubRow(_.Y,X.h.i)
     });
     WriteLn("TRACE "+_MID+"5");
-    Set X.weight = get.X.weight(X.dif);
-    Set X.weight.der = get.X.weight.partial.derivative.1(X.dif,s);
+    Set X.weight = get.X.weight(X.dif, X.h_);
+    Set X.weight.der = get.X.weight.partial.derivative.1(X.dif,X.h_,s);
     WriteLn("TRACE "+_MID+"6");
     VMatrix YnCum = DifEq(1/(1-B),_.Yn);