restart; with(plots):
 P := [ [0,0], [1,4], [3,5], [4,0]]; # Control points for a Bezier curve
 f := factor(expand((1-t)^3*P[1] + 3*(1-t)^2*t*P[2] + 3*(1-t)*t^2*P[3] + t^3*P[4]));
 fp := diff(f,t);
 fpp := diff(fp,t);
 ds := sqrt(factor(fp[1]^2+fp[2]^2));
 T := factor(expand(fp/ds));  #Unit tangent vector
 Ns := [-T[2],T[1]]; # Signed normal
 ks := (-fpp[1]*fp[2] + fpp[2]*fp[1])/ds^3; #Signed curvature
 C := factor(expand(f+(1/ks)*Ns)); #Center of circle of curvature
 R := abs(1/ks); #Radius of curvature
 a := 0; b := 1;  n := 32;
 P := plot([f[1],f[2], t=a..b],color=green,thickness=2): #Plot the original
 for i from 0 to n do  tt := a + (b-a)*i/n;  xx := evalf( subs(t=tt,C)); RR := evalf(subs(t=tt,R)); 
   Tmp := plot([xx[1]+RR*cos(theta),xx[2]+RR*sin(theta),theta=0..2*Pi],color=blue,numpoints=100): #Plot of osculating circle at tt
   Q[i] := display([P,Tmp]);
 od:

 display([Q[j] $j=0..n],insequence=true,scaling=constrained,view=[-5..5,-5..5]);