@mmcdara Interesting and useful information. Thank you!

I've seen even library procedures using:

if cond or is(cond) then ... fi

@vv We can generate infinitely many solutions e.g. playing with u1 and u3 and using fsolve.

solfs:=fsolve([perim, u1=8, u3=12]); # infinitely many solutions, playing e.g. with u1 and u3
   #    {u1 = 8.000000000, u2 = 7.152963125, u3 = 12.00000000,  v1 = 0.5346258468, v2 = 0.2316337764, v3 = 0.5061472158,  v4 = 0.1464801560, v5 = 0.3671814714, v6 = 0.4712668716}

eval(seq("A"||i=A||i,i=0..8), solfs);
#   "A0" = [0, 0], "A1" = [8.000000000, 0],    "A2" = [7.152963125, 12.00000000], "A3" = [4.277006774, 0],   "A4" = [7.803797650, 2.779605317], 
#   "A5" = [7.571274644, 6.073766590],  "A6" = [1.047767154, 1.757761872],  "A7" = [2.626435525, 4.406177657],  "A8" = [4.231660510, 2.239322836]

eval([perim], solfs);
#    [10.00000000 = 10, 15.00000000 = 15, 11.00000000 = 11,
#      9.000000001 = 9, 13.00000000 = 13, 12.00000000 = 12, 
#      20.00000000 = 20]

plots:-display(plot(eval([A0,A1,A2,A0],solfs)), plot(eval([A3,A4,A6,A3],solfs)), plot(eval([A5,A7,A8,A5],solfs)), scaling=constrained);

So, you want to prove that if X in M_n(R) commutes with any invertible matrix A then X is a scalar multiple of Id (= identity matrix).

(Here R is a unitary and commutative ring).

This is easy to prove by taking A = Id + E, where E has a zero diagonal and a single nonzero entry. 
But not so easy with Maple, and I don't understand why  you want it; it would be very artificial and useless.

@one man  Here is a compact version for animation:

restart;

with(plots):
A := [cos(t),sin(t)]:
sol := solve({x^2+y^2-1=0,  (x-cos(t))^2+(y-sin(t))^2 - sin(t)^2 =0}, [x,y], explicit):
P1:=rhs~(sol[1]): P2:=rhs~(sol[2]):
F := proc(a, dt:=0.02) Threads:-Sleep(dt);  display(
  plot([A[],t=0..2*Pi]), 
  plot(eval([P1,P2], t=a), symbolsize=20, color=red),
  plot(eval([A, [cos(t),0]], t=a), symbolsize=20, color=gold, thickness=4),
  plot([eval([P1,P2], t=a)], color=red, thickness=4), 
  plot([cos(t),sin(t)/2, t=0..a], color=blue, thickness=4),  
  implicitplot( (x-cos(a))^2+(y-sin(a))^2 - sin(a)^2 =0, x=-2..2,y=-2..2, color=green),
  textplot([[cos(a),sin(a), "A"], [cos(a),sin(a)/2, "C"]],align=left),
  scaling=constrained, view=[-2..2, -2..2]   ) end proc:

Explore( F(a), a=0..2*Pi, animate, autorun, loop );

The fact that the locus of C is an ellipse is easy to prove:

restart;
A := [cos(t),sin(t)]:
sol := solve({x^2+y^2-1=0,  (x-cos(t))^2+(y-sin(t))^2 - sin(t)^2 =0}, [x,y], explicit):
P1:=rhs~(sol[1]): P2:=rhs~(sol[2]):
C:=simplify(P1+lambda*~(P2-P1)):
lambda := solve(C[1]=cos(t), lambda):
C:=simplify(C);

              C := [cos(t), (1/2)*sin(t)]

So, the locus is an ellipse with semiaxes 1 and 1/2.

@mmcdara I never recommend to students to solve such problems using just Maple. But it is nice to know how to use Maple to check maths resuts while learning.

@Carl Love Yes, you are right

eps:=1e-5:
Optimization:-LPSolve(AD, {CD=AB-2, BC=AB+2, AD=AB+BC+CD, AB>=eps, BC>=eps, CD>=eps}, integervariables={AD});
#   [7, [AB = 2.33333333333333, AD = 7, BC = 4.33333333333333,  CD = 0.333333333333333]]

Or, using exact ("pure") maths:

inf(eval(AD::integer, solve({CD=AB-2, BC=AB+2, AD=AB+BC+CD, AB=a, BC=b, CD=c})), a>0,b>0,c>0);

#  inf((6 + 3*c)::integer, 0 < a, 0 < b, 0 < c)  ;   ==>  7

@JAMET Your code is for conics with center. A parabola has no center.

@one man Yes, it's about precision.
It is easy to exhibit a circle with r = 0.7698144869481772563348519  (Digits=25) having a unique tangency point. But visually is already obvious.

sol := [0.75, [r = 0.75, x0 = 1.08, x1 = 0.699729170630851, x2 = 1.27446720511737, x3 = 1.89574141571275, y0 = 1.16, y1 = 0.505162888466810, y2 = 1.89574141571275, y3 = 1.27446720511737]]:
p1:=plots:-implicitplot(e, 0..2, 0..2): p2:=plots:-implicitplot(h, 0..2, 0..2):
p3:=plottools:-circle(eval([x0,y0],sol[2]), eval(r,sol[2]), color=red):
p4:=plots:-pointplot(eval([[x0,y0],[x1,y1],[x2,y2],[x3,y3]],sol[2]), symbolsize=8, color=blue):
p5:=plots:-textplot(eval({seq}([x||i,y||i,P__||i], i=0..3), sol[2]), align={below, left}, color=blue):
plots:-display(p1,p2,p3,p4,p5);

(I modified sol manually. I have not computed, but visually it seems I am right.)

@one man  If you impose a unique point of tangency on each curve, then there is no circle with maximal radius. 
In fact, for any r' < r = 0.7698... (and r' > 0.75, say) we can decrease the radius of the circle to r'  and move it a bit to be tangent just near P₁ and near P₂ .

@nm I don't think so. You simply apply it outside your context.

Let's say you have a function F defined on sets, but F({g,{h}})  is not necessarily equal to F({g,h}).
Then, if you apply Carl's FlattenSet to

r:={a,{b,c},d,{e,f,F({g,{h}})}};

you will come up here too with "This code has serious bug in it".

This is a problem in number theory, not quite a Diophantine equation; an IMO level one, so not very simple.
It can be transformed into a Diophantine equation:

isolve( x^2+y^2 = k*(1+x*y));

          {k = _Z1^2, x = _Z1, y = 0}

Actually this is given by Maple 2018 (on an old Laptop of mine). This is not correct without x>=0, y>=0,  for example x=2, y=-1, k=-5.
In Maple 2024 the result is NULL (also wrong, of course).

Note that isolve is not even able to solve directly the equation for k=4,  although it is of a Pell type. But we can help here:

isolve(z^2-3*y^2 - 4 = 0);  # x = 2*y + z

  ==>  z= ..., y=...

So, unfortunately isolve is not useful for our problem. The only thing we can do is to check numerically the result:

for i to 1000 do for j to i-1 do
  if irem(  ( i^2+j^2) , 1+i*j)=0 then lprint([i,j],(i^2+j^2)/(1+i*j)) fi
od od:

[8, 2], 4
[27, 3], 9
[30, 8], 4
[64, 4], 16
[112, 30], 4
[125, 5], 25
[216, 6], 36
[240, 27], 9
[343, 7], 49
[418, 112], 4
[512, 8], 64
[729, 9], 81
[1000, 10], 100

@mmcdara You have chosen the HermiteH polynomials in the Galerkin method. They are orthogonal in the interval (-oo, oo) with the weight exp(-x^2) you did not use.
So, why HermiteH? Probably ChebyshevT with the corresponding weight are better.