read bezier:
with (plots):

BezierApprox := proc(x,t)
	local i,n,p;
	n := nops(x);
	p := sum(x[i]*B(n-1,i-1,t), i=1..n);
end:

PlotPoly := proc(Px,Py,x,y,T)
	local plot1, plot2, px, py;
	px := subs(X=x,Px);
	py := subs(Y=y,Py);
	plot1 := plot([px,py, t=0..1]):
	plot2 := plot([cx,cy,t=0..1]):
	print(plot({px, cx}, t=0..1, title=`Approximation to x component of circle`));
	print(plot({py, cy}, t=0..1, title=`Approximation to y component of circle`));
	display({plot1,plot2},title=T);
end:

SubsPoly := proc(px,py,x,y)
	plot([subs(X=x,px), subs(Y=y,py), t=0..1]);
end:

x := 0,X,1,1;
y := 1,1,Y,0;

px := BezierApprox([x],t):
py := BezierApprox([y],t):
cx := 2*t/(1+t^2):
cy := (1-t^2)/(1+t^2):

dx := px - cx;
dy := py - cy;

ix := int( dx^2, t=0..1);
iy := int( dy^2, t=0..1);
ry := Re(solve(iy=0,Y)[1]);
rx := Re(solve(ix=0,X)[1]);
print(`P2 is `, (evalf(rx),1));
print(`P3 is `, (1,evalf(ry)));
PlotPoly(px,py, rx,ry, `Optimal Cubic Bezier Approximation to Circle`);

p1 := SubsPoly(px,py, 1.0, 0.0):
p2 := SubsPoly(px,py, 0.9, 0.1):
p3 := SubsPoly(px,py, 0.8, 0.2):
p4 := SubsPoly(px,py, 0.7, 0.3):
p5 := SubsPoly(px,py, 0.6, 0.4):
p6 := SubsPoly(px,py, 0.5, 0.5):
p7 := SubsPoly(px,py,  rx,  ry):
p8 := plot([cx,cy, t=0..1], colour=GREEN):
display({p1,p2,p3,p4,p5,p6,p7,p8},title=`Some Random Cubic Bezier Approximations to Circle`);