> restart: with(plots):
> decas:=proc(degree, coeff,t)
> # de Casteljau-Algorithmus zur Berechnung einer Koordinate der
> Bezierkurve
> # Input: degree: Grad der Kurve
> #        coeff: Feld der Kontrollpunkte (entsprechende Koordinate)
> #        t: Parameter
> # option trace;
> local r,i, t1, coeffa;
> t1:=1-t; 
> for i from 1 to degree+1 do coeffa[i]:=coeff[i] od;
> for r from 1 to degree do
>     for i from 0 to degree -r do coeffa[i+1]:=t1*coeffa[i+1]+
> t*coeffa[i+2] od
> od;
> evalf(coeffa[1])
> end:
> bezier:=proc(degree, points, it)
> # Plotten einer parametrischen Bezierkurve mit Kontrollpunkten points
> # it : Anzahl der Iterationen zur Teilung des Kontrollpolygons
> # option trace;
> global Curve;
> local i,k,r,coneu,coalt,coneu1, coeff, bezier_curve, cpoints, subdiv;
> subdiv:=proc(degree,coord, t)
> # Berechnung des linken und rechten Teil-Kontrollpolygons fuer t mit
> de Casteljau
> # Abspeicherung der neuen Kontrollpunkte in coneu
> # option trace;
> local coleft, coright, t1, i,r, coneu;
> t1:=1-t; coright:=vector(degree+1); coleft:=vector(degree+1);
> coneu:=vector(2*degree+1);
> for i from 0 to degree  do  coright[i+1]:=coord[i+1] od;
> for r from 1 to degree do
>    for i from 0 to degree - r do
>       coright[i+1]:=t1*coright[i+1] + t*coright[i+2]
>    od
> od;
> for i from 0 to degree  do  coleft[degree-i+1]:=coord[i+1] od;
> for r from 1 to degree do
>    for i from 0 to degree - r do
>       coleft[i+1]:=t*coleft[i+1] + t1*coleft[i+2]
>    od
> od;
> for i from 1 to degree+1 do coneu[i]:=coleft[degree+2-i];
> coneu[i+degree]:=coright[i] od;
> evalm(coneu)
> end:
> # x Komponente
> Curve:=array(1..2^(it)*degree+1,1..2);
> for i from 1 to degree +1 do coneu[i]:=points[i,1] od;
> for i from 1 to it do
>   for r from 1 to  2^(i-1)*degree +1 do coalt[r]:=coneu[r] od;
>   for k from 0 to 2^(i-1)-1 do 
>      for r from 1 to degree +1 do coeff[r]:=coalt[r+k*degree] od;
>      coneu1:= subdiv(degree,coeff,0.5);
>      for r from 1 to 2*degree +1 do coneu[r+ 2*k*degree]:=coneu1[r]
> od;
>   od;
> od;
> for i from 1 to 2^(it)*degree +1 do Curve[i,1]:=coneu[i] od;
> # y-Komponente
> for i from 1 to degree +1 do coneu[i]:=points[i,2] od;
> for i from 1 to it do
>   for r from 1 to  2^(i-1)*degree +1 do coalt[r]:=coneu[r] od;
>   for k from 0 to 2^(i-1)-1 do 
>      for r from 1 to degree +1 do coeff[r]:=coalt[r+k*degree] od;
>      coneu1:= subdiv(degree,coeff,0.5);
>      for r from 1 to 2*degree +1 do coneu[r+ 2*k*degree]:=coneu1[r]
> od;
>   od;
> od;
> for i from 1 to 2^(it)*degree+1 do Curve[i,2]:=coneu[i] od;
> bezier_curve:=plot(convert(Curve,listlist)):
> cpoints:=pointplot(points,color=blue,symbol=circle):
> RETURN (display([bezier_curve,cpoints])):
> end:
> P:=[[1,2],[2,6],[4,4],[5,3]];
> display([bezier(3,P,0), bezier(3,P,1), bezier(3,P,4)]);

                P := [[1, 2], [2, 6], [4, 4], [5, 3]]


> P:=[[1,1],[3,3],[2,-5],[-2,-3],[-1,3]];

          P := [[1, 1], [3, 3], [2, -5], [-2, -3], [-1, 3]]

> display([bezier(4,P,0),bezier(4,P,1),bezier(4,P,3)]);

>