The inscribed square problem, also known as the Toeplitz conjecture, is an unsolved quastion in geometry: Does every plane simple closed curve (Jordan curve) contain all four vertices of some square? This is true if the curve is convex or piecewise smooth and in other special cases. The problem was proposed by Otto Toeplitz in 1911. For detailes see  https://en.wikipedia.org/wiki/Inscribed_square_problem

The Inscribed_Square procedure finds numerically one or more solutions for a curve defined by parametric equations of its boundary or by the equation F(x,y)=0. The required parameter of procedure  L  is the list of equations of the boundary links or the equation  F(x,y)=0 . Optional parameters:  N  and  R . By default  N='onesolution' (the procedure finds one solution), if  N  is any symbol (for example  N='s'), then more solutions.  R  is the range for the length of the side of the square (by defalt  R=0.1..100 ).

The second procedure  Pic  visualizes the results obtained.

The codes of the procedures:

restart;
Inscribed_Square:=proc(L::{list(list),`=`},N::symbol:='onesolution',R::range:=0.1..100)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, M, eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, sol, Sol;
uses LinearAlgebra;
if L::list then
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
M:=<0,-1;1,0>;
eq1:=eval(C[1])=eval((B+M.(B-A))[1]);
eq2:=eval(C[2])=eval((B+M.(B-A))[2]);
eq3:=eval(D[1])=eval((C+M.(C-B))[1]);
eq4:=eval(D[2])=eval((C+M.(C-B))[2]);
eq5:=eval(DotProduct(B-A,B-A, conjugate=false))=d^2;
sol:=fsolve([eq1,eq2,eq3,eq4,eq5],{t1=op([2,2,1],L0[i])..op([2,2,2],L0[i]),t2=op([2,2,1],L0[j])..op([2,2,2],L0[j]),t3=op([2,2,1],L0[k])..op([2,2,2],L0[k]),t4=op([2,2,1],L0[m])..op([2,2,2],L0[m]),d=R});
if type(sol,set(`=`)) then if N='onesolution' then return convert~(eval([A,B,C,D],sol),list) else c:=c+1; Sol[c]:=convert~(eval([A,B,C,D],sol),list) fi;
 fi; 
od: od: od: od:
Sol:=fnormal(convert(Sol,list),7);
print(Sol);
ListTools:-Categorize((X,Y)->`and`(seq(is(convert(X,set)[i]=convert(Y,set)[i]),i=1..4)) , Sol);
return map(t->t[1],[%]);
else
A,B,C,D:=<x1,y1>,<x2,y2>,<x3,y3>,<x4,y4>:
M:=<0,-1;1,0>:
eq1:=eval(C[1])=eval((B+M.(B-A))[1]):
eq2:=eval(C[2])=eval((B+M.(B-A))[2]):
eq3:=eval(D[1])=eval((C+M.(C-B))[1]):
eq4:=eval(D[2])=eval((C+M.(C-B))[2]):
eq5:=eval(LinearAlgebra:-DotProduct((B-A,B-A), conjugate=false))=d^2:
eq6:=eval(L,[x=x1,y=y1]):
eq7:=eval(L,[x=x2,y=y2]):
eq8:=eval(L,[x=x3,y=y3]):
eq9:=eval(L,[x=x4,y=y4]):
sol:=fsolve({eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9},{seq([x||i=-2..2,y||i=-2..2][],i=1..4),d=R});
eval([[x1,y1],[x2,y2],[x3,y3],[x4,y4]], sol):
fi;
end proc:

Pic:=proc(L,Sol,R::range:=-20..20)
local P1, P2, P3, T;
uses plots, plottools;
P1:=`if`(L::list,seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L), implicitplot(L, x=R,y=R, color=blue, thickness=2, gridrefine=3));
P2:=polygon(Sol,color=yellow,thickness=0);
P3:=curve([Sol[],Sol[1]],color=red,thickness=3):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
display(P1,P2,P3,T, scaling=constrained, size=[800,500], axes=none);
end proc:

Examples of use:

The curve consists of a semicircle, a segment and a semi-ellipse (find 1 solution):

L:=[[[cos(t),sin(t)],t=0..Pi],[[t,0],t=-1..0],[[0.5+0.5*cos(t),0.8*sin(t)],t=Pi..2*Pi]]:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

The procedure finds 6 solutions for a non-convex pentagon:

 L:=[[[0,0],[9,0]],[[9,0],[8,5]],[[8,5],[5,3]],[[5,3],[0,4]],[[0,4],[0,0]]]:
Sol:=Inscribed_Square(L,'s');
plots:-display(Matrix(3,2,[seq(Pic(L,Sol[i]),i=1..6)]),size=[300,200]);

For an implicitly defined curve, only one solution can be found:

L:=abs(x)+2*abs(y)-sin((2*x-y))-cos(x+y)^2=3:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

               
See more examples in the attached file.

Inscribed_Square.mw