> 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]);