restart; with(plots):

Warning, the name changecoords has been redefined

P := [ [0,0], [1,4], [3,5], [4,0]]; # Control points for a Bezier curve

NiM+SSJQRzYiNyY3JCIiIUYoNyQiIiIiIiU3JCIiJCIiJjckRitGKA==

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

NiM+SSJmRzYiNyQsJComSSJ0R0YlIiIiLCgqJEYpIiIjRi1GKSEiJEYuRipGKiEiIiwkKihGKUYqLCZGKUYqIiIlRipGKiwmRilGKkYvRipGKkYu

fp := diff(f,t);

NiM+SSNmcEc2IjckLCoqJEkidEdGJSIiIyEiI0YpIiIkRiwiIiIqJkYpRi0sJkYpIiIlISIkRi1GLSEiIiwoKiYsJkYpRi1GMEYtRi0sJkYpRi1GMkYtRi1GMSomRilGLUY2Ri1GMSomRilGLUY1Ri1GMQ==

fpp := diff(fp,t);

NiM+SSRmcHBHNiI3JCwmSSJ0R0YlISM3IiInIiIiLCZGKCEjPUYtRis=

ds := sqrt(factor(fp[1]^2+fp[2]^2));

NiM+SSNkc0c2IiwkKiQsLCokSSJ0R0YlIiIlIiM4KiRGKiIiJCIjR0YqISNXIiM8IiIiKiRGKiIiIyIjNyNGMkY0Ri4=

T := factor(expand(fp/ds)); #Unit tangent vector

NiM+SSJURzYiNyQsJComLCgqJEkidEdGJSIiI0YsRishIiMhIiIiIiJGLywsKiRGKyIiJSIjOCokRisiIiQiI0dGKyEjVyIjPEYvRioiIzcjRi5GLEYuLCQqJiwoRipGNUYrIiInISIlRi9GL0YwRjpGLg==

Ns := [-T[2],T[1]]; # Signed normal

NiM+SSNOc0c2IjckKiYsKCokSSJ0R0YlIiIjIiIkRioiIichIiUiIiJGLywsKiRGKiIiJSIjOCokRipGLCIjR0YqISNXIiM8Ri9GKSIjNyMhIiJGKywkKiYsKEYpRitGKiEiI0Y6Ri9GL0YwRjlGOg==

ks := (-fpp[1]*fp[2] + fpp[2]*fp[1])/ds^3; #Signed curvature

NiM+SSNrc0c2IiwkKiYsJiomLCZJInRHRiUhIzciIiciIiJGLiwoKiYsJkYrRi4iIiVGLkYuLCZGK0YuISIiRi5GLiEiJComRitGLkYzRi5GNSomRitGLkYxRi5GNUYuRjQqJiwmRishIz1GOkYuRi4sKiokRisiIiMhIiNGKyIiJEY/Ri4qJkYrRi4sJkYrRjJGNUYuRi5GNEYuRi5GLiwsKiRGK0YyIiM4KiRGK0Y/IiNHRishI1ciIzxGLkY8IiM3I0Y1Rj0jRi4iI0Y=

C := factor(expand(f+(1/ks)*Ns)); #Center of circle of curvature

NiM+SSJDRzYiNyQsJComLDAqJEkidEdGJSIiJyIkPCIqJEYrIiImIiRBJiokRisiIiQhJDcmKiRGKyIiIyEkJnoqJEYrIiIlIiQjUUYrIiQjeiEkLyMiIiJGPCwoRjQiIipGKyEiJiIiKEY8ISIiI0ZBRjUsJComLDBGKiIjRUYuIiM3RjEhI19GNCIjR0Y3ISMqKUYrIiNtISM8RjxGPEY9RkEjRjJGNQ==

R := abs(1/ks); #Radius of curvature

NiM+SSJSRzYiLCQtSSRhYnNHSSpwcm90ZWN0ZWRHRik2IyomLCYqJiwmSSJ0R0YlISM3IiInIiIiRjIsKComLCZGL0YyIiIlRjJGMiwmRi9GMiEiIkYyRjIhIiQqJkYvRjJGN0YyRjkqJkYvRjJGNUYyRjlGMkY4KiYsJkYvISM9Rj5GMkYyLCoqJEYvIiIjISIjRi8iIiRGQ0YyKiZGL0YyLCZGL0Y2RjlGMkYyRjhGMkYyRjgsLCokRi9GNiIjOCokRi9GQyIjR0YvISNXIiM8RjJGQCIjNyNGQ0ZBIiNG

a := 0; b := 1; n := 32; P := plot([f[1],f[2], t=a..b],color=green,thickness=2): #Plot the originalfor 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:

NiM+SSJhRzYiIiIhNiM+SSJiRzYiIiIiNiM+SSJuRzYiIiNL

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

Error, invalid plot structure