restat;
with(plots);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
a := 11;
b := 5;
theta := Pi/5;

ell := (x*cos(theta) + y*sin(theta))^2/a^2 + (x*sin(theta) - y*cos(theta))^2/b^2 = 1;
sol := solve({ell, y = m}, {x, y}, explicit);
subs(sol[1], x) = subs(sol[2], x);
sol := solve(%, m);
bb := sol[1];
evalf(%);
                          -7.626741393

sol := solve({ell, x = m}, {x, y}, explicit);
subs(sol[1], y) = subs(sol[2], y);
sol := solve(%, m);

aa := sol[1];
evalf(%);
                          -9.371916332

A := [aa, bb];
B := [aa, -bb];
C := [-aa, -bb];
Dd := [-aa, bb];
Ell := implicitplot(ell, x = -12 .. 12, y = -12 .. 12, thickness = 3, color = blue);
quadri := plot([A, B, C, Dd, A], color = pink, filled = true);

display([quadri, Ell], scaling = constrained, axes = none);

As "it is this program works;I wish it rotates when theta varies.

Very good program; Thanks Tom Leslie .And without  using Geometry ?

restart:local D:
  with(plots):
  with(plottools):
 

a := 7;
b := a/2;
NULL;

r := b/2;
c := 2*sqrt(10)*r;
A := [c/2, c/2];
B := [c/2, -1/2*c];
C := [-c/2, (-c)/2];
D := [-c/2, c/2];
f := (x, y) -> x^2/a^2 + y^2/b^2 = 1;

Ell1 := implicitplot(5/98*x^2 - 3/49*x*y + 5/98*y^2 = 1, x = -9 .. 9, y = -6 .. 6, color = white);
g := (x, y) -> x^2/b^2 + y^2/a^2 = 1;

Ell2 := implicitplot(5/98*x^2 + 3/49*x*y + 5/98*y^2 = 1, x = -9 .. 9, y = -6 .. 6, color = white);
quadri := plot([A, B, C, D, A], x = -6 .. 6, y = -6 .. 6, filled = true, color = red, transparency = 0.6);
Cir1 := disk([c/2 - r, -c/2 + r], r, color = blue);
Cir2 := disk([-c/2 + r, -c/2 + r], r, color = blue);
Cir3 := disk([-c/2 + r, c/2 - r], r, color = blue);
Cir4 := disk([c/2 - r, c/2 - r], r, color = blue);
display([quadri, Ell1, Ell2, Cir1, Cir2, Cir3, Cir4], scaling = constrained, axes = normal);

Is it possible to have while the inside of the ellipses ?

Thanks a lot.
Here is my last pgr.
How to animate this drawing moving P1 in the circle dfor instance

restart;
with(plots): with(plottools):with(LinearAlgebra): unprotect(D);
_EnvHorizontalName := 'x':
_EnvVerticalName := 'y':
Vdot := proc (U, V)local i: add(U[i]*V[i], i = 1 .. 2) end proc:
dist := proc (M, N) sqrt(Vdot(expand(M-N), expand(M-N))) end proc:
MinDistPoint := proc (A::[algebraic, algebraic], B::[algebraic, algebraic], P::[algebraic, algebraic]) 
local D, R, V; 
description "Point on line AB at minimum distance from P";
D := `<,>`(A-B); R := `<,>`(`<|>`(0, -1), `<|>`(1, 0)); 
V := `<,>`(`<,>`(P) . D, -R . `<,>`(A) . `<,>`(B)); `~`[`/`]([V . D, R . V . D], ` $`, D . D) 
end proc:

EqBIS := proc(P, U, V) local a, eq1, M1, t, PU, PV, bissec1; a := (P - U)/LinearAlgebra:-Norm(P - U, 2) + (P - V)/LinearAlgebra:-Norm(P - V, 2); M1 := P + a*t; eq1 := op(eliminate({x = M1[1], y = M1[2]}, t)); RETURN(op(eq1[2])); end proc;

Cen := proc(M, N, R) local eq1, eq2, sol; eq1 := EqBIS(M, N, R) = 0; eq2 := EqBIS(N, M, R) = 0; sol := simplify(solve({eq1, eq2}, {x, y})); RETURN([subs(sol, x), subs(sol, y)]); end proc;

CircleParm := t -> [(-t^2 + 1)/(t^2 + 1), 2*t/(t^2 + 1)];
NULL;
P1 := Vector(CircleParm(1/3)):

P2 := Vector(CircleParm(5)):
P3 := Vector(CircleParm(-1/10)):
P4 := Vector(CircleParm(-19/2)):
p1 := convert(P1, list);
p2 := convert(P2, list);
p3 := convert(P3, list);
p4 := convert(P4, list);
C4 := Vector(Cen(P1, P2, P3)):H := convert(C4, list):
C1 := Vector(Cen(P4, P2, P3));
J := convert(C1, list);
C2 := Vector(Cen(P4, P1, P3));
Ii := convert(C2, list);
C3 := Vector(Cen(P4, P2, P1));
K := convert(C3, list);
Pts := [P1, P2, P3, P4, C1, C2, C3, C4]:
r4 := dist(MinDistPoint(p1, p2, H), H);
r2 := dist(MinDistPoint(p1, p4, Ii), Ii);
r1 := dist(MinDistPoint(p2, p3, J), J);
r3 := dist(MinDistPoint(p2, p1, K), K);
Cir1 := circle(Ii, r2, color = black);
Cir2 := circle(K, r3, color = black);
Cir3 := circle(J, r1, color = black);
Cir4 := circle(H, r4, color = black);
Poly := polygonplot([p1, p2, p4, p3, p1], color = blue, transparency = 0.9);
SegP1P4 := plot([p1, p4], color = black);
SegP2P3 := plot([p2, p3], color = black);
quadri := plot([p1, p2, p4, p3, p1], color = black);
rect := polygonplot([H, Ii, J, K, H], color = red, transparency = 0.7);
tex := textplot([[p1[], "P1"], [p2[], "P2"], [p3[], "P3"], [p4[], "P4"], [H[], "H"], [Ii[], "I"], [J[], "J"], [K[], "K"]], align = ["above", "right"]);
display(Poly, tex, quadri, rect, SegP1P4, SegP2P3, Cir1, Cir2, Cir3, Cir4, implicitplot([x^2 + y^2 - 1], x = -2 .. 2, y = -4 .. 2, colour = [blue], scaling = constrained), pointplot(Pts, symbolsize = 16), axes = none);

how to draw all the lines P1P2,P2P3,P3P4, etc,  and the incirles of centers C1, C2..

Thanks a lot, Still an error. 
Tri1 := plot([P2, P3, P4, P2], color = red);
Error, (in plot) incorrect first argument [Vector(2, {(1) = -12/13, (2) = 5/13}), Vector(2, {(1) = 99/101, (2) = -20/101}), Vector(2, {(1) = -5/13, (2) = -12/13}), Vector(2, {(1) = -12/13, (2) = 5/13})]

Sorry.
I have still error message :
  Error, (in plots:-display) expecting plot structures but received: [plots:-pointplot*([seq*([cos((2/17)*k*Pi), sin((2/17)*k*Pi)], k = 1 .. 16)], symbol = solidcircle, color = red, symbolsize = 10), plots:-textplot*([seq*([cos((2/17)*k*Pi), sin((2/17)*k*Pi), "Mk"], k = 1 .. 16)], align = ["above", "right"])]

How to use textplot to write M1...M17 near each point ? Thank you.

clearly how to show that the only Friday 13 of 2022 is in May. Thank you very much.

Thank you a lot for your answer.
One observation : good explanations are given in internet site : WolframMathWorld Pell equation.

pellsolve := proc(D::posint)
local P, Q, a, A, B, i; if type(sqrt(D), integer) then error "D must be a nonsquare integer"; end if;
P := 0; Q := 1; a := floor(sqrt(D)); A := 1, a; B := 0, 1;
for i do P := a*Q - P;
Q := (D - P^2)/Q;
a := floor((P + sqrt(D))/Q);
A := A[2], a*A[2] + A[1]; B := B[2], a*B[2] + B[1]; break; if Q = 1 and i mod 2 = 0;
end do; return A[1], B[1]; end proc;
pellsolve(28);for example
####retourne la solution minimale (x,y) de l'équation de Pell x`^2 - Dy = 1, Fonctionne bien.

I’m surprised I can’t find the same coordinates using the barycentric coordinates !

Vdot := proc (U, V)
local i: add(U[i]*V[i], i = 1 .. 2) end:
dist := proc (M, N) sqrt(Vdot(expand(M-N), expand(M-N))) end:
NULL;
A := [0, 0]:
B := [7, 0]:
C := [3, 5]:
a := dist(B, C):
b := dist(A, C):
c := dist(A, B):
An := arccos((-a^2 + b^2 + c^2)/(2*b*c));
Bn := arccos((a^2 - b^2 + c^2)/(2*a*c));
Cn := arccos((a^2 + b^2 - c^2)/(2*a*b));
evalf(An + Bn + Cn);
(A[1]*csc(An + Pi/3) + B[1]*csc(Bn + Pi/3) + C[1]*csc(Cn + Pi/3))/(csc(An + Pi/3) + csc(Bn + Pi/3) + csc(Cn + Pi/3)), (A[2]*csc(An + Pi/3) + B[2]*csc(Bn + Pi/3) + C[2]*csc(Cn + Pi/3))/(csc(An + Pi/3) + csc(Bn + Pi/3) + csc(Cn + Pi/3));
evalf(%);
Do you see why ? Thank you.

Thanks a lot. Now how to locate vertex and focus of the parabola in the base (O,i,j) ?

Very good for the order of the coefficients but there is still a confusion between parabola and ellipse. How to avoid it ? Thank you.

Thank you very much. Without Geometry.