diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 7b662acb7084f1042e09c35a119e288db1f5ded0..732f87a34eb2f4516156b301a058791836010575 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -135,6 +135,32 @@ SPoint2 InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
   }
   return p;
 }
+// Bezier
+static Vertex InterpolateBezier(Curve *Curve, double u, int derivee)
+{
+  int NbCurves = (List_Nbr(Curve->Control_Points) - 1) / 3;
+  int iCurve = (int)floor(u * (double)NbCurves);
+  if(iCurve == NbCurves) iCurve -= 1; // u = 1
+  double t1 = (double)(iCurve) / (double)(NbCurves);
+  double t2 = (double)(iCurve+1) / (double)(NbCurves);
+  double t = (u - t1) / (t2 - t1);
+  Vertex *v[4];
+  for(int i = 0; i < 4; i++) {
+    List_Read(Curve->Control_Points, iCurve * 3 + i , &v[i]);
+  }
+
+  if(Curve->geometry){
+    SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, t1, t2, Curve->geometry,derivee);
+    SPoint3 pt = Curve->geometry->point(pp);
+    Vertex V;
+    V.Pos.X = pt.x();
+    V.Pos.Y = pt.y();
+    V.Pos.Z = pt.z();
+    return V;
+  }
+  else
+    return InterpolateCubicSpline(v, t, Curve->mat, derivee, t1, t2);
+}
 
 // Uniform BSplines
 static Vertex InterpolateUBS(Curve *Curve, double u, int derivee)
@@ -274,10 +300,13 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
   if(derivee==2) {
     switch (c->Typ) {
     case MSH_SEGM_BSPLN:
-    case MSH_SEGM_BEZIER:
       V = InterpolateUBS(c, u, 2);
       V.u = u;
       break;
+    case MSH_SEGM_BEZIER:
+      V = InterpolateBezier(c, u, 2);
+      V.u = u;
+      break;
     default :
       double eps1 = (u == 0) ? 0 : 1.e-5;
       double eps2 = (u == 1) ? 0 : 1.e-5;
@@ -358,9 +387,11 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
     break;
 
   case MSH_SEGM_BSPLN:
-  case MSH_SEGM_BEZIER:
     V = InterpolateUBS(c, u, 0);
     break;
+  case MSH_SEGM_BEZIER:
+    V = InterpolateBezier(c, u, 0);
+    break;
 
   case MSH_SEGM_NURBS:
     V = InterpolateNurbs(c, u, 0);
diff --git a/demos/splines.geo b/demos/splines.geo
index 27b86c2f83290f051c2657fdcfac7112c288238f..53b2cc650b7caba67e8e248286c0ea5c9f66220f 100644
--- a/demos/splines.geo
+++ b/demos/splines.geo
@@ -56,8 +56,7 @@ Point(p+1) = {xx+.5, 0,  0, lc} ;
 Point(p+2) = {xx+.5, 1, 0, 0.1*lc} ;
 Point(p+3) = {xx, 1, 0, lc} ;
 l = newreg;
-//Bezier(l) = {p+3,p+2,p+1,p}; // Bezier curves are broken
-BSpline(l) = {p+3,p+2,p+1,p};
+Bezier(l) = {p+3,p+2,p+1,p}; // Bezier curves are broken
 Line(l+1) = {p,p+3};
 s = newreg;
 Line Loop(s) = {-l,-(l+1)};