diff --git a/dox/user_guides/draw_test_harness/draw_test_harness.md b/dox/user_guides/draw_test_harness/draw_test_harness.md
index dc31956d97..d4be9c003e 100644
--- a/dox/user_guides/draw_test_harness/draw_test_harness.md
+++ b/dox/user_guides/draw_test_harness/draw_test_harness.md
@@ -5552,17 +5552,20 @@ intersect e c p
 
 Syntax:      
 ~~~~~
-2dintersect curve1 curve2 
+2dintersect curve1 [curve2] [-tol tol] [-state]
 ~~~~~
 
-Displays the intersection points between two 2d curves. 
+Displays the intersection points between 2d curves.
+Options:
+ -tol - allows changing the intersection tolerance (default value is 1.e-3);
+ -state - allows printing the intersection state for each point.
 
 **Example:** 
 ~~~~~
 # intersect two 2d ellipses 
 ellipse e1 0 0 5 2 
 ellipse e2 0 0 0 1 5 2 
-2dintersect e1 e2 
+2dintersect e1 e2 -tol 1.e-10 -state
 ~~~~~
 
 @subsubsection occt_draw_6_7_3 intconcon
diff --git a/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx b/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx
index ff430f6115..40e2b417de 100644
--- a/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx
+++ b/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx
@@ -645,6 +645,9 @@ void BOPAlgo_PaveFiller::AnalyzeShrunkData(const Handle(BOPDS_PaveBlock)& thePB,
         AddWarning (new BOPAlgo_AlertTooSmallEdge (aWarnShape));
       else
         AddWarning (new BOPAlgo_AlertBadPositioning (aWarnShape));
+      Standard_Real aTS1, aTS2;
+      theSR.ShrunkRange(aTS1, aTS2);
+      thePB->SetShrunkData(aTS1, aTS2, Bnd_Box(), Standard_False);
       return;
     }
     //
diff --git a/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx b/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx
index 930e1a82d8..36e1c9fa72 100644
--- a/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx
+++ b/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx
@@ -1474,6 +1474,58 @@ Standard_Boolean BOPAlgo_PaveFiller::IsExistingPaveBlock
   }
   return bRet;
 }
+
+//=======================================================================
+//function : 
+//purpose  : 
+//=======================================================================
+static void getBoundPaves(const BOPDS_DS* theDS,
+                          BOPDS_Curve& theNC,
+                          Standard_Integer theNV[2])
+{
+  theNV[0] = theNV[1] = -1;
+
+  // get extreme paves
+  Handle(BOPDS_PaveBlock)& aPB = theNC.ChangePaveBlock1();
+  const BOPDS_ListOfPave& aLP = aPB->ExtPaves();
+  Standard_Integer aNbEP = aLP.Extent();
+  if (aNbEP == 0)
+    return;
+  Standard_Real aTmin = RealLast();
+  Standard_Real aTmax = -aTmin;
+  for (BOPDS_ListIteratorOfListOfPave aItLP(aLP); aItLP.More(); aItLP.Next())
+  {
+    const BOPDS_Pave& aPv = aItLP.Value();
+    Standard_Integer nV;
+    Standard_Real aTV;
+    aPv.Contents(nV, aTV);
+    if (aTV < aTmin) {
+      theNV[0] = aPv.Index();
+      aTmin = aTV;
+    }
+    if (aTV > aTmax) {
+      theNV[1] = aPv.Index();
+      aTmax = aTV;
+    }
+  }
+
+  // compare extreme vertices with ends of the curve
+  const IntTools_Curve& aIC = theNC.Curve();
+  Standard_Real aT[2];
+  gp_Pnt aP[2];
+  aIC.Bounds(aT[0], aT[1], aP[0], aP[1]);
+  Standard_Real aTol = Max(theNC.Tolerance(), theNC.TangentialTolerance());
+  aTol += Precision::Confusion();
+  for (Standard_Integer j = 0; j < 2; ++j)
+  {
+    const BOPDS_ShapeInfo& aSIV = theDS->ShapeInfo(theNV[j]);
+    const TopoDS_Vertex& aV = (*(TopoDS_Vertex *)(&aSIV.Shape()));
+    Standard_Integer iFlag = BOPTools_AlgoTools::ComputeVV(aV, aP[j], aTol);
+    if (iFlag != 0)
+      theNV[j] = -1;
+  }
+}
+
 //=======================================================================
 //function : PutBoundPaveOnCurve
 //purpose  : 
@@ -1483,95 +1535,50 @@ void BOPAlgo_PaveFiller::PutBoundPaveOnCurve(const TopoDS_Face& aF1,
                                              BOPDS_Curve& aNC,
                                              TColStd_ListOfInteger& aLVB)
 {
-  Standard_Boolean bVF;
-  Standard_Integer nV, iFlag, nVn, j, aNbEP;
-  Standard_Real aT[2], aTmin, aTmax, aTV, aTol, aTolVnew;
-  gp_Pnt aP[2];
-  TopoDS_Vertex aVn;
-  BOPDS_ListIteratorOfListOfPave aItLP;
-  BOPDS_Pave aPn, aPMM[2];
-  //
-  aTolVnew = Precision::Confusion();
-  //
   const IntTools_Curve& aIC=aNC.Curve();
+  Standard_Real aT[2];
+  gp_Pnt aP[2];
   aIC.Bounds(aT[0], aT[1], aP[0], aP[1]);
   Standard_Real aTolR3D = Max(aNC.Tolerance(), aNC.TangentialTolerance());
+  Handle(BOPDS_PaveBlock)& aPB = aNC.ChangePaveBlock1();
+  // Get numbers of vertices assigned to the ends of the curve
+  Standard_Integer aBndNV[2];
+  getBoundPaves(myDS, aNC, aBndNV);
   //
-  Handle(BOPDS_PaveBlock)& aPB=aNC.ChangePaveBlock1();
-  const BOPDS_ListOfPave& aLP=aPB->ExtPaves();
-  //
-  aNbEP=aLP.Extent();
-  if (aNbEP) {
-    aTmin=1.e10;
-    aTmax=-aTmin;
-    //
-    aItLP.Initialize(aLP);
-    for (; aItLP.More(); aItLP.Next()) {
-      const BOPDS_Pave& aPv=aItLP.Value();
-      aPv.Contents(nV, aTV);
-      if (aTV<aTmin) {
-        aPMM[0]=aPv;
-        aTmin=aTV;
+  Standard_Real aTolVnew = Precision::Confusion();
+  for (Standard_Integer j = 0; j<2; ++j)
+  {
+    if (aBndNV[j] < 0)
+    {
+      // no vertex on this end
+      if (j && aP[1].IsEqual(aP[0], aTolVnew)) {
+        //if curve is closed, process only one bound
+        continue;
       }
-      if (aTV>aTmax) {
-        aPMM[1]=aPv;
-        aTmax=aTV;
-      }
-    }
-  }
-  //
-  for (j=0; j<2; ++j) {
-    //if curve is closed, process only one bound
-    if (j && aP[1].IsEqual(aP[0], aTolVnew)) {
-      continue;
-    }
-    //
-    iFlag=1;
-    //
-    if (aNbEP) {
-      Bnd_Box aBoxP;
-      //
-      aBoxP.Set(aP[j]);
-      aTol = aTolR3D+Precision::Confusion();
-      aBoxP.Enlarge(aTol);
-      const BOPDS_Pave& aPV=aPMM[j];
-      nV=aPV.Index();
-      const BOPDS_ShapeInfo& aSIV=myDS->ShapeInfo(nV);
-      const TopoDS_Vertex& aV=(*(TopoDS_Vertex *)(&aSIV.Shape()));
-      const Bnd_Box& aBoxV=aSIV.Box();
-      if (!aBoxP.IsOut(aBoxV)){
-        iFlag=BOPTools_AlgoTools::ComputeVV(aV, aP[j], aTol);
-      }
-    }
-    if (iFlag) {
-      // 900/L5
-      bVF=myContext->IsValidPointForFaces (aP[j], aF1, aF2, aTolR3D);
+      Standard_Boolean bVF = myContext->IsValidPointForFaces(aP[j], aF1, aF2, aTolR3D);
       if (!bVF) {
         continue;
       }
-      //
-      BOPDS_ShapeInfo aSIVn;
-      //
+      TopoDS_Vertex aVn;
       BOPTools_AlgoTools::MakeNewVertex(aP[j], aTolR3D, aVn);
+      BOPTools_AlgoTools::UpdateVertex(aIC, aT[j], aVn);
+      aTolVnew = BRep_Tool::Tolerance(aVn);
+
+      BOPDS_ShapeInfo aSIVn;
       aSIVn.SetShapeType(TopAbs_VERTEX);
       aSIVn.SetShape(aVn);
-      //
-      nVn=myDS->Append(aSIVn);
-      //
+
+      Bnd_Box& aBox = aSIVn.ChangeBox();
+      BRepBndLib::Add(aVn, aBox);
+      aBox.SetGap(aBox.GetGap() + Precision::Confusion());
+
+      Standard_Integer nVn = myDS->Append(aSIVn);
+
+      BOPDS_Pave aPn;
       aPn.SetIndex(nVn);
       aPn.SetParameter(aT[j]);
       aPB->AppendExtPave(aPn);
-      //
-      aVn=(*(TopoDS_Vertex *)(&myDS->Shape(nVn)));
-      BOPTools_AlgoTools::UpdateVertex (aIC, aT[j], aVn);
-      //
-      aTolVnew = BRep_Tool::Tolerance(aVn);
-      //
-      BOPDS_ShapeInfo& aSIDS=myDS->ChangeShapeInfo(nVn);
-      Bnd_Box& aBoxDS=aSIDS.ChangeBox();
-      BRepBndLib::Add(aVn, aBoxDS);
-      aBoxDS.SetGap(aBoxDS.GetGap() + Precision::Confusion());
-      //
+
       aLVB.Append(nVn);
     }
   }
@@ -1913,7 +1920,7 @@ void BOPAlgo_PaveFiller::GetEFPnts
 }
 
 //=======================================================================
-//function : ProcessUnUsedVertices
+//function : PutStickPavesOnCurve
 //purpose  : 
 //=======================================================================
   void BOPAlgo_PaveFiller::PutStickPavesOnCurve
@@ -1925,6 +1932,14 @@ void BOPAlgo_PaveFiller::GetEFPnts
    TColStd_DataMapOfIntegerReal& aMVTol,
    TColStd_DataMapOfIntegerListOfInteger& aDMVLV)
 {
+  // Get numbers of vertices assigned to the ends of the curve
+  Standard_Integer aBndNV[2];
+  getBoundPaves(myDS, aNC, aBndNV);
+  if (aBndNV[0] >= 0 && aBndNV[1] >= 0)
+  {
+    // both curve ends already have assigned vertices
+    return;
+  }
   TColStd_MapOfInteger aMV;
   aMV.Assign(aMVStick);
   RemoveUsedVertices(aNC, aMV);
@@ -1961,6 +1976,8 @@ void BOPAlgo_PaveFiller::GetEFPnts
       aPV=BRep_Tool::Pnt(aV);
       //
       for (m=0; m<2; ++m) {
+        if (aBndNV[m] >= 0)
+          continue;
         aD2=aPC[m].SquareDistance(aPV);
         if (aD2>aDT2) {// no rich
           continue; 
@@ -1990,8 +2007,6 @@ void BOPAlgo_PaveFiller::GetEFPnts
       }
     }//for (jVU=1; jVU=aNbVU; ++jVU) {
   }
-  //}//if (aTypeC==GeomAbs_BezierCurve || aTypeC==GeomAbs_BSplineCurve) {
-  //}//if(aType1==GeomAbs_Torus  || aType2==GeomAbs_Torus) {
 }
 
 //=======================================================================
@@ -2706,11 +2721,9 @@ void BOPAlgo_PaveFiller::UpdatePaveBlocks
         if (wasRegularEdge && !isDegEdge && nV[0] == nV[1]) {
           // now edge has the same vertex on both ends;
           // check if it is not a regular closed curve.
-          const TopoDS_Edge& aE = TopoDS::Edge(myDS->Shape(nE));
-          const TopoDS_Vertex& aV = TopoDS::Vertex(myDS->Shape(nV[0]));
-          Standard_Real aLength = IntTools::Length(aE);
-          Standard_Real aTolV = BRep_Tool::Tolerance(aV);
-          if (aLength <= aTolV * 2.) {
+          FillShrunkData(aPB);
+          if (!aPB->HasShrunkData())
+          {
             // micro edge, so mark it for removal
             aMicroEdges.Add(nE);
             continue;
diff --git a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx
index 11f4daa960..89e5f30bea 100644
--- a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx
+++ b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx
@@ -301,53 +301,89 @@ static Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer n, const
 //=======================================================================
 static Standard_Integer intersect(Draw_Interpretor& di, Standard_Integer n, const char** a)
 {
-  if( n < 2) 
+  if (n < 2)
   {
-    cout<< "2dintersect curve curve [Tol]"<<endl;
+    di.PrintHelp(a[0]);
     return 1;
   }
 
-  Standard_Integer k = 1;
-  Handle(Geom2d_Curve) C1 = DrawTrSurf::GetCurve2d(a[k++]);
-  if ( C1.IsNull()) 
+  Handle(Geom2d_Curve) C1 = DrawTrSurf::GetCurve2d(a[1]);
+  if (C1.IsNull())
+  {
+    di << "Curve " << a[1] << " is null\n";
     return 1;
-
-  Standard_Real Tol = 0.001;
-  Geom2dAPI_InterCurveCurve Intersector;
+  }
 
   Handle(Geom2d_Curve) C2;
-  if ( k < n ) {
-    C2 = DrawTrSurf::GetCurve2d(a[k++]);
-    if ( C2.IsNull())
-      return 1;
-  }
-  if(k < n)
-    Tol = Draw::Atof(a[k]);
+  Standard_Real Tol = 0.001;
+  Standard_Boolean bPrintState = Standard_False;
 
-  if(!C2.IsNull())
+  // Retrieve other parameters if any
+  for (Standard_Integer i = 2; i < n; ++i)
   {
-    Intersector.Init(C1,C2,Tol);
+    if (!strcmp(a[i], "-tol"))
+    {
+      Tol = Draw::Atof(a[++i]);
+    }
+    else if (!strcmp(a[i], "-state"))
+    {
+      bPrintState = Standard_True;
+    }
+    else
+    {
+      C2 = DrawTrSurf::GetCurve2d(a[i]);
+      if (C2.IsNull())
+      {
+        di << "Curve " << a[i] << " is null\n";
+        return 1;
+      }
+    }
   }
-  else {
+
+  Geom2dAPI_InterCurveCurve Intersector;
+
+  if (!C2.IsNull())
+    // Curves intersection
+    Intersector.Init(C1, C2, Tol);
+  else
+    // Self-intersection of the curve
     Intersector.Init(C1, Tol);
+
+  const Geom2dInt_GInter& anIntTool = Intersector.Intersector();
+  if (!anIntTool.IsDone())
+  {
+    di << "Intersection failed\n";
+    return 0;
   }
 
-  Standard_Integer i;
+  if (anIntTool.IsEmpty())
+    return 0;
 
-  for ( i = 1; i <= Intersector.NbPoints(); i++) {
+  Standard_Integer aNbPoints = Intersector.NbPoints();
+  for (Standard_Integer i = 1; i <= aNbPoints; i++)
+  {
+    // API simplified result
     gp_Pnt2d P = Intersector.Point(i);
-
-    di<<"Intersection point "<<i<<" : "<<P.X()<<" "<<P.Y()<<"\n";
-    di<<"parameter on the fist: "<<Intersector.Intersector().Point(i).ParamOnFirst();
-    di<<" parameter on the second: "<<Intersector.Intersector().Point(i).ParamOnSecond()<<"\n";
-    Handle(Draw_Marker2D) mark = new Draw_Marker2D( P, Draw_X, Draw_vert); 
+    di << "Intersection point " << i << " : " << P.X() << " " << P.Y() << "\n";
+    // Intersection extended results from intersection tool
+    const IntRes2d_IntersectionPoint& aPInt = anIntTool.Point(i);
+    di << "parameter on the fist: " << aPInt.ParamOnFirst();
+    di << " parameter on the second: " << aPInt.ParamOnSecond() << "\n";
+    if (bPrintState)
+    {
+      di << "Intersection type: " <<
+        (aPInt.TransitionOfFirst().IsTangent() ? "TOUCH" : "INTERSECTION") << "\n";
+    }
+    Handle(Draw_Marker2D) mark = new Draw_Marker2D(P, Draw_X, Draw_vert);
     dout << mark;
   }
   dout.Flush();
 
-  Handle(Geom2d_Curve) S1,S2;
+  Handle(Geom2d_Curve) S1, S2;
   Handle(DrawTrSurf_Curve2d) CD;
-  for ( i = 1; i <= Intersector.NbSegments(); i++) {
+  Standard_Integer aNbSegments = Intersector.NbSegments();
+  for (Standard_Integer i = 1; i <= aNbSegments; i++)
+  {
     di << "Segment #" << i << " found.\n";
     Intersector.Segment(i,S1,S2);
     CD = new DrawTrSurf_Curve2d(S1, Draw_bleu, 30);
@@ -540,8 +576,13 @@ void GeomliteTest::API2dCommands(Draw_Interpretor& theCommands)
 
   g = "GEOMETRY intersections";
 
-  theCommands.Add("2dintersect", "intersect curve curve [Tol]",__FILE__,
-		  intersect,g);
+  theCommands.Add("2dintersect", "2dintersect curve1 [curve2] [-tol tol] [-state]\n"
+                   "Intersects the given 2d curve(s)."
+                   "If only one curve is given, it will be checked on self-intersection.\n"
+                   "Options:\n"
+                   " -tol - allows changing the intersection tolerance (default value is 1.e-3);\n"
+                   " -state - allows printing the intersection state for each point.",
+                   __FILE__, intersect, g);
 
   theCommands.Add("2dintanalytical", "intersect circle1 and circle2 using IntAna",__FILE__,
 		  intersect_ana,g);
diff --git a/src/math/math_FunctionRoots.cxx b/src/math/math_FunctionRoots.cxx
index 7e19758b51..a6bb20c694 100644
--- a/src/math/math_FunctionRoots.cxx
+++ b/src/math/math_FunctionRoots.cxx
@@ -22,6 +22,7 @@
 #include <math_DirectPolynomialRoots.hxx>
 #include <math_FunctionRoots.hxx>
 #include <math_FunctionWithDerivative.hxx>
+#include <math_BracketedRoot.hxx>
 #include <Standard_RangeError.hxx>
 #include <StdFail_NotDone.hxx>
 #include <TColStd_Array1OfReal.hxx>
@@ -37,6 +38,21 @@ static Standard_Boolean myDebug = 0;
 static Standard_Integer nbsolve = 0;
 #endif
 
+class DerivFunction: public math_Function
+{
+  math_FunctionWithDerivative *myF;
+
+public:
+  DerivFunction(math_FunctionWithDerivative& theF)
+    : myF (&theF)
+  {}
+
+  virtual Standard_Boolean Value(const Standard_Real theX, Standard_Real& theFval)
+  {
+    return myF->Derivative(theX, theFval);
+  }
+};
+
 static void  AppendRoot(TColStd_SequenceOfReal& Sol,
 			TColStd_SequenceOfInteger& NbStateSol,
 			const Standard_Real X,
@@ -415,78 +431,126 @@ math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
 	  }
 	}
 	if(Rediscr) {
-	  //-- --------------------------------------------------
-	  //-- On recherche un extrema entre x0 et x3
-	  //-- x1 et x2 sont tels que x0<x1<x2<x3 
-	  //-- et |f(x0)| > |f(x1)|   et |f(x3)| > |f(x2)|
-	  //--
-	  //-- En entree : a=xm-dx  b=xm c=xm+dx
-	  Standard_Real x0,x1,x2,x3,f0,f3;
-	  Standard_Real R=0.61803399;
-	  Standard_Real C=1.0-R;
-	  Standard_Real tolCR=NEpsX*10.0;
-	  f0=ptrval(im1);
-	  f3=ptrval(ip1);
-	  x0=xm-dx;
-	  x3=xm+dx;
-          if(x0 < X0) x0=X0;
-          if(x3 > XN) x3=XN;
-	  Standard_Boolean recherche_minimum = (f0>0.0);
+          Standard_Real x0 = xm - dx;
+          Standard_Real x3 = xm + dx;
+          if (x0 < X0) x0 = X0;
+          if (x3 > XN) x3 = XN;
+          Standard_Real aSolX1 = 0., aSolX2 = 0.;
+          Standard_Real aVal1 = 0., aVal2 = 0.;
+          Standard_Real aDer1 = 0., aDer2 = 0.;
+          Standard_Boolean isSol1 = Standard_False;
+          Standard_Boolean isSol2 = Standard_False;
+          //-- ----------------------------------------------------
+          //-- Find minimum of the function |F| between x0 and x3
+          //-- by searching for the zero of the function derivative
+          DerivFunction aDerF(F);
+          math_BracketedRoot aBR(aDerF, x0, x3, EpsX);
+          if (aBR.IsDone())
+          {
+            aSolX1 = aBR.Root();
+            F.Value(aSolX1, aVal1);
+            aVal1 = Abs(aVal1);
+            if (aVal1 < EpsF)
+            {
+              isSol1 = Standard_True;
+              aDer1 = aBR.Value();
+            }
+          }
 
-	  if(Abs(x3-xm) > Abs(x0-xm)) { x1=xm; x2=xm+C*(x3-xm); } 
-	  else                        { x2=xm; x1=xm-C*(xm-x0); } 
-	  Standard_Real f1,f2;
-	  F.Value(x1,f1); f1-=K;
-	  F.Value(x2,f2); f2-=K;
-	  //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
-	  Standard_Real tolX = 0.001 * NEpsX;
-	  while(Abs(x3-x0) > tolCR*(Abs(x1)+Abs(x2)) && (Abs(x1 -x2) > tolX)) { 
-	    //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ", 
-	    //--    x0,f0,x1,f1,x2,f2,x3,f3);
-	    if(recherche_minimum) {  
-	      if(f2<f1) { 
-		x0=x1; x1=x2; x2=R*x1+C*x3;
-		f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
-	      }
-	      else { 
-		x3=x2; x2=x1; x1=R*x2+C*x0;
-		f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
-	      }
-	    }
-	    else { 
-	      if(f2>f1) { 
-		x0=x1; x1=x2; x2=R*x1+C*x3;
-		f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
-	      }
-	      else { 
-		x3=x2; x2=x1; x1=R*x2+C*x0;
-		f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
-	      }
-	    }
-	    //-- On ne fait pas que chercher des extremas. Il faut verifier 
-	    //-- si on ne tombe pas sur une racine 
-	    if(f1*f0 <0.0) { 
-	      //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
-	      Solve(F,K,x0,f0,x1,f1,tol,NEpsX,Sol,NbStateSol);
-	    }
-	    if(f2*f3 <0.0)  { 
-	      //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
-	      Solve(F,K,x2,f2,x3,f3,tol,NEpsX,Sol,NbStateSol);
-	    }
-	  }
-	  if(f1<f2) { 
-	    //-- x1,f(x1) minimum
-	    if(Abs(f1) < EpsF) { 
-	      AppendRoot(Sol,NbStateSol,x1,F,K,NEpsX);
-	    }
-	  }
-	  else { 
-	    //-- x2.f(x2) minimum
-	    if(Abs(f2) < EpsF) { 
-	      AppendRoot(Sol,NbStateSol,x2,F,K,NEpsX);
-	    }
-	  }
-	} //-- Recherche d un extrema    
+          //-- --------------------------------------------------
+          //-- On recherche un extrema entre x0 et x3
+          //-- x1 et x2 sont tels que x0<x1<x2<x3 
+          //-- et |f(x0)| > |f(x1)|   et |f(x3)| > |f(x2)|
+          //--
+          //-- En entree : a=xm-dx  b=xm c=xm+dx
+          Standard_Real x1, x2, f0, f3;
+          Standard_Real R = 0.61803399;
+          Standard_Real C = 1.0 - R;
+          Standard_Real tolCR = NEpsX*10.0;
+          f0 = ptrval(im1);
+          f3 = ptrval(ip1);
+          Standard_Boolean recherche_minimum = (f0 > 0.0);
+
+          if (Abs(x3 - xm) > Abs(x0 - xm)) { x1 = xm; x2 = xm + C*(x3 - xm); }
+          else                        { x2 = xm; x1 = xm - C*(xm - x0); }
+          Standard_Real f1, f2;
+          F.Value(x1, f1); f1 -= K;
+          F.Value(x2, f2); f2 -= K;
+          //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
+          Standard_Real tolX = 0.001 * NEpsX;
+          while (Abs(x3 - x0) > tolCR*(Abs(x1) + Abs(x2)) && (Abs(x1 - x2) > tolX)) {
+            //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ", 
+            //--    x0,f0,x1,f1,x2,f2,x3,f3);
+            if (recherche_minimum) {
+              if (f2 < f1) {
+                x0 = x1; x1 = x2; x2 = R*x1 + C*x3;
+                f0 = f1; f1 = f2; F.Value(x2, f2); f2 -= K;
+              }
+              else {
+                x3 = x2; x2 = x1; x1 = R*x2 + C*x0;
+                f3 = f2; f2 = f1; F.Value(x1, f1); f1 -= K;
+              }
+            }
+            else {
+              if (f2 > f1) {
+                x0 = x1; x1 = x2; x2 = R*x1 + C*x3;
+                f0 = f1; f1 = f2; F.Value(x2, f2); f2 -= K;
+              }
+              else {
+                x3 = x2; x2 = x1; x1 = R*x2 + C*x0;
+                f3 = f2; f2 = f1; F.Value(x1, f1); f1 -= K;
+              }
+            }
+            //-- On ne fait pas que chercher des extremas. Il faut verifier 
+            //-- si on ne tombe pas sur une racine 
+            if (f1*f0 < 0.0) {
+              //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
+              Solve(F, K, x0, f0, x1, f1, tol, NEpsX, Sol, NbStateSol);
+            }
+            if (f2*f3 < 0.0)  {
+              //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
+              Solve(F, K, x2, f2, x3, f3, tol, NEpsX, Sol, NbStateSol);
+            }
+          }
+          if ((recherche_minimum && f1<f2) || (!recherche_minimum && f1>f2)) {
+            //-- x1,f(x1) minimum
+            if (Abs(f1) < EpsF) {
+              isSol2 = Standard_True;
+              aSolX2 = x1;
+              aVal2 = Abs(f1);
+            }
+          }
+          else {
+            //-- x2.f(x2) minimum
+            if (Abs(f2) < EpsF) {
+              isSol2 = Standard_True;
+              aSolX2 = x2;
+              aVal2 = Abs(f2);
+            }
+          }
+          // Choose the best solution between aSolX1, aSolX2
+          if (isSol1 && isSol2)
+          {
+            if (aVal2 - aVal1 > EpsF)
+              AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
+            else if (aVal1 - aVal2 > EpsF)
+              AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
+            else
+            {
+              aDer1 = Abs(aDer1);
+              F.Derivative(aSolX2, aDer2);
+              aDer2 = Abs(aDer2);
+              if (aDer1 < aDer2)
+                AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
+              else
+                AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
+            }
+          }
+          else if (isSol1)
+            AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
+          else if(isSol2)
+            AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
+        } //-- Recherche d un extrema    
       } //-- for     
     }      
     
diff --git a/tests/bugs/modalg_5/bug23706_43 b/tests/bugs/modalg_5/bug23706_43
index 1751d4615e..a4d29de49f 100644
--- a/tests/bugs/modalg_5/bug23706_43
+++ b/tests/bugs/modalg_5/bug23706_43
@@ -11,15 +11,13 @@ set y 8.0
 set z -2.0
 set pp_ch1 1.0371228345434986
 set pp_ch2 0.99999999682789309
-set pp_ch3 0.99999999486742297
 
 restore [locate_data_file bug23706_c07.draw] c
 set info [proj c $x $y $z]
 
 regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
 regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
-regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
-if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
    puts "Error : Projection is not correct"
 } else {
    puts "OK: Projection is correct"
diff --git a/tests/bugs/modalg_5/bug23706_48 b/tests/bugs/modalg_5/bug23706_48
index 45d7e629d6..4302bae90a 100644
--- a/tests/bugs/modalg_5/bug23706_48
+++ b/tests/bugs/modalg_5/bug23706_48
@@ -11,15 +11,13 @@ set y 8.0
 set z -2.0
 set pp_ch1 1.0371228345434986
 set pp_ch2 0.99999999851019361
-set pp_ch3 1.0000000000000002
 
 restore [locate_data_file bug23706_c08.draw] c
 set info [proj c $x $y $z]
 
 regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
 regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
-regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
-if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
    puts "Error : Projection is not correct"
 } else {
    puts "OK: Projection is correct"
diff --git a/tests/bugs/modalg_6/bug25593 b/tests/bugs/modalg_6/bug25593
index 40db85cd9d..9917a0d5c4 100755
--- a/tests/bugs/modalg_6/bug25593
+++ b/tests/bugs/modalg_6/bug25593
@@ -10,12 +10,12 @@ restore [locate_data_file bug25593_face1.brep] b
 pcurve b
 
 puts "Intersection 1"
-set info1 [2dintersect b_2 b_3 1.e-10]
+set info1 [2dintersect b_2 b_3 -tol 1.e-10]
 regexp {Intersection point 1 : +([-0-9.+eE]+) +([-0-9.+eE]+)} $info1 full p1x p1y
 regexp {parameter on the fist: +([-0-9.+eE]+) +parameter on the second: +([-0-9.+eE]+)} $info1 full par1f par1s
 
 puts "Intersection 2"
-set info2 [2dintersect b_3 b_2 1.e-10]
+set info2 [2dintersect b_3 b_2 -tol 1.e-10]
 regexp {Intersection point 1 : +([-0-9.+eE]+) +([-0-9.+eE]+)} $info2 full p2x p2y
 regexp {parameter on the fist: +([-0-9.+eE]+) +parameter on the second: +([-0-9.+eE]+)} $info2 full par2f par2s
 
diff --git a/tests/bugs/modalg_6/bug26938_3 b/tests/bugs/modalg_6/bug26938_3
index 42ced6af77..ba99301605 100644
--- a/tests/bugs/modalg_6/bug26938_3
+++ b/tests/bugs/modalg_6/bug26938_3
@@ -1,4 +1,3 @@
-puts "TODO OCC24694 All: Error : The result of cut operation is self-interfered shape"
 puts "TODO OCC24694 ALL: Result is WRONG because number of SHELL entities in shape \"result\" is 2"
 puts "TODO OCC24694 ALL: Result is WRONG because number of SOLID entities in shape \"result\" is 2"
 puts "=========="
diff --git a/tests/bugs/modalg_6/bug26938_4 b/tests/bugs/modalg_6/bug26938_4
index 2ac00095dc..44c0575fdd 100644
--- a/tests/bugs/modalg_6/bug26938_4
+++ b/tests/bugs/modalg_6/bug26938_4
@@ -1,4 +1,3 @@
-puts "TODO OCC24694 ALL: Error : The result of cut operation is self-interfered shape"
 puts "TODO OCC24694 ALL: Result is WRONG because number of SHELL entities in shape \"result\" is 2"
 puts "TODO OCC24694 ALL: Result is WRONG because number of SOLID entities in shape \"result\" is 2"
 
diff --git a/tests/bugs/modalg_7/bug28274 b/tests/bugs/modalg_7/bug28274
index 752b172c1d..5581707de2 100755
--- a/tests/bugs/modalg_7/bug28274
+++ b/tests/bugs/modalg_7/bug28274
@@ -1,5 +1,3 @@
-puts "TODO OCC28274 ALL: Error: BOPAlgo_MakerVolume fails to create solid"
-
 puts "========"
 puts "OCC28274"
 puts "========"
@@ -14,12 +12,12 @@ bnondestructive 1
 
 set Log [eval mkvolume result [explode a f] ]
 
-if {[regexp "Warning" ${Log}] == 1} {
-    puts "Error: BOPAlgo_MakerVolume fails to create solid"
-}
-
 checkshape result
 
-bopcheck result
+if ![regexp "OK" [bopcheck result]] {
+  puts "Error: result has self-interferred shapes"
+}
 
+checknbshapes result -solid 3 -shell 3 -face 11 -wire 11
+checkprops result -s 985.692 -v 741.787
 checkview -display result -2d -path ${imagedir}/${test_image}.png
diff --git a/tests/bugs/modalg_7/bug28385_11 b/tests/bugs/modalg_7/bug28385_11
new file mode 100644
index 0000000000..8cdc8c4b5e
--- /dev/null
+++ b/tests/bugs/modalg_7/bug28385_11
@@ -0,0 +1,17 @@
+puts "========"
+puts "OCC28385"
+puts "========"
+puts ""
+#################################################
+# Improve drawing isolines (DBRep_IsoBuilder algorithm)
+#################################################
+
+
+binrestore [locate_data_file bug28385_de_iges_1_J9_faces.bin] a
+
+smallview
+isos a 2
+repeat 9 u
+fit
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png
diff --git a/tests/bugs/moddata_2/bug23587 b/tests/bugs/moddata_2/bug23587
index 2e083838eb..ee41d7906d 100755
--- a/tests/bugs/moddata_2/bug23587
+++ b/tests/bugs/moddata_2/bug23587
@@ -17,10 +17,10 @@ pcurve c1 a_1 f1
 pcurve c2 a_2 f1
 pcurve c3 a_3 f1
 pcurve c4 a_4 f1
-set inter1 [2dintersect c1 c2 1.e-6]
-set inter2 [2dintersect c3 c2 1.e-6]
-set inter3 [2dintersect c1 c4 1.e-6]
-set inter4 [2dintersect c3 c4 1.e-6]
+set inter1 [2dintersect c1 c2 -tol 1.e-6]
+set inter2 [2dintersect c3 c2 -tol 1.e-6]
+set inter3 [2dintersect c1 c4 -tol 1.e-6]
+set inter4 [2dintersect c3 c4 -tol 1.e-6]
 set int1 [regexp {Intersection point 1} $inter1]
 set int2 [regexp {Intersection point 1} $inter2]
 set int3 [regexp {Intersection point 1} $inter3]
diff --git a/tests/bugs/moddata_3/bug25876 b/tests/bugs/moddata_3/bug25876
index b4e1062e15..826485789c 100755
--- a/tests/bugs/moddata_3/bug25876
+++ b/tests/bugs/moddata_3/bug25876
@@ -15,7 +15,7 @@ explode e2
 mk2dcurve c1 e1_1
 mk2dcurve c2 e2_1
 
-set inter1 [2dintersect c1 c2 1e-3]
+set inter1 [2dintersect c1 c2 -tol 1e-3]
 set int1 [regexp {Intersection point 2} ${inter1}]
 if { ${int1} == 0 } {
     puts "Error : Second intersection is not found"
@@ -27,7 +27,7 @@ view 1 -2D- 728 450 400 400
 xwd ${imagedir}/${test_image}_1.png
 donly c1 c2
 
-set inter2 [2dintersect c1 c2 1e-7]
+set inter2 [2dintersect c1 c2 -tol 1e-7]
 set int2 [regexp {Intersection point 2} ${inter2}]
 if { ${int2} == 0 } {
     puts "Error : Second intersection is not found"
diff --git a/tests/lowalgos/2dinter/A1 b/tests/lowalgos/2dinter/A1
new file mode 100644
index 0000000000..c77cf96749
--- /dev/null
+++ b/tests/lowalgos/2dinter/A1
@@ -0,0 +1,8 @@
+restore [locate_data_file bug29368_1.brep] a
+explode a
+mk2dcurve c1 a_1 1
+mk2dcurve c2 a_2 1
+set out [2dintersect c1 c2 -tol 1.e-9 -state]
+if ![regexp "TOUCH" $out] {
+  puts "Error: intersection type is not TOUCH"
+}
diff --git a/tests/lowalgos/2dinter/A2 b/tests/lowalgos/2dinter/A2
new file mode 100644
index 0000000000..439dbd627d
--- /dev/null
+++ b/tests/lowalgos/2dinter/A2
@@ -0,0 +1,8 @@
+restore [locate_data_file bug29368_2.brep] a
+explode a
+mk2dcurve c1 a_1 1
+mk2dcurve c2 a_2 1
+set out [2dintersect c1 c2 -tol 1.e-9 -state]
+if ![regexp "TOUCH" $out] {
+  puts "Error: intersection type is not TOUCH"
+}
diff --git a/tests/lowalgos/begin b/tests/lowalgos/begin
new file mode 100644
index 0000000000..b752ad2bbc
--- /dev/null
+++ b/tests/lowalgos/begin
@@ -0,0 +1,14 @@
+if { [array get Draw_Groups "TOPOLOGY Check commands"] == "" } {
+    pload TOPTEST
+    pload AISV
+}
+# To prevent loops limit to 1 minutes
+cpulimit 60
+
+if { [info exists imagedir] == 0 } {
+   set imagedir .
+}
+
+if { [info exists test_image ] == 0 } {
+   set test_image photo
+}
diff --git a/tests/lowalgos/end b/tests/lowalgos/end
new file mode 100644
index 0000000000..57901e2abe
--- /dev/null
+++ b/tests/lowalgos/end
@@ -0,0 +1,2 @@
+# to end a test script
+puts "TEST COMPLETED"
diff --git a/tests/lowalgos/grids.list b/tests/lowalgos/grids.list
new file mode 100644
index 0000000000..8b614e18bc
--- /dev/null
+++ b/tests/lowalgos/grids.list
@@ -0,0 +1 @@
+001 2dinter
diff --git a/tests/lowalgos/parse.rules b/tests/lowalgos/parse.rules
new file mode 100644
index 0000000000..610d206403
--- /dev/null
+++ b/tests/lowalgos/parse.rules
@@ -0,0 +1 @@
+FAILED /\bFaulty\b/ bad shape