Actually, in Maple it is possible to define the vertices of the polygon simply by clicking! You have to use the so called Embedded components.

For an example, just type ?examples,ProgrammaticContentGeneration and execute the commands in this worksheet to see how it works (press Enter several times and click a few times in the plot when it appears).

@acer  As I remember, years ago, sometimes, it could be necessary to install some version of 
Visual C++ Redistributable for Visual Studio.

@Alfred_F  I use 1d input and 2d output. Maple promotes now 2d input; it looks nicer but may introduce hidden characters hard to detect, so, non-working code. 2d input has some limitations too; note that the 2d notation is internally converted to 1d. I prefer the old school.

@Alfred_F You are using 2d input; just replace the last line with:

lprint(cat("That's all for x <=y <= z <= ", N))

are quickly found by Maple:

restart;
s1:=x+y+z: s2:=x*y+x*z+y*z: s3:=x*y*z:
N:=1000:
for x to N do for y from x to N do for z from y to N do
  if isqrt(s1)^2=s1 and isqrt(s2)^2=s2 and isqrt(s3)^2=s3
     then lprint([x,y,z]) fi;
od od od:
lprint("That's all for x <=y <= z <= " || N);

[45, 64, 180]
[65, 104, 360]
[72, 136, 153]
[80, 225, 320]
[98, 272, 306]
[180, 256, 720]
[272, 306, 578]
[288, 544, 612]
That's all for x <=y <= z <= 1000

An interesting fact mentioned by Broman (see my answer) is that Holditch's curve may not exist if the length is not short enough. An example is the ellipse with semiaxes 3 and 1 and the chord length 2, considered by @dharr. In the article:
Monterde J.,  Rochera D. - Holditch's Ellipse Unveiled. The American Mathematical Monthly , Vol. 124, No. 5 (May 2017), pp. 403-421,
a maximal  chord length is computed for an ellipse with semiaxes p,q: 
 

@Alfred_F The beauty of the result you mentioned in the title is generated by the theoretical fact, not by some animations for a chord on a curve (easy to obtain).

restart
Digits:=15:
X:=t -> 10*cos(t) + 7*exp(sin(t))/3:
Y:=t -> 7*sin(t)-sin(sin(t))/2:
a, b, r, n := 4, 3, 3, 100:
tB:=proc(tA)
  local t;
  fsolve( (X(tA)-X(t))^2 + (Y(tA)-Y(t))^2 = (a+b)^2, t=tA..tA+r);
end proc:
vt:=evalf([seq](2*Pi/n*k,k=0..n)):
vA :=[seq]([X(t),Y(t)],t=vt):
vB := [seq]( [ X,Y](tB(t)), t=vt):
vC := a/(a+b)*~vA + b/(a+b)*~vB:
area:= proc(v) local n:=nops(v), i;
 (add(v[i,1]*v[i+1,2] - v[i+1,1]*v[i,2], i=1..n-1) + v[n,1]*v[1,2] - v[1,1]*v[n,2])/2;
end proc:
da:=(area(vA) - area(vC))/(a*b): 
with(plots):
pAC := display(plot(vA), plot(vC), scaling=constrained, caption=typeset("area/(a*b)=",da)):
animate(plot, [ [vA[floor(i)], vB[floor(i)]] ,color=gold, thickness=4], i=1..n+1, background=pAC,  scaling=constrained);

Edit: added the area computation.

@Alfred_F @Alfred_F  Yes, but then the proof is not completely rigurous.
Here is the original proof (1854) via:
Broman A. - Holditch's Theorem, Mathematics Magazine, Vol. 54, No. 3 (May, 1981), pp. 99-108.

@sand15 What you say is valid in the French strict formalism (Bourbaki style, which I like very much, probably because I am not young :-)). But the French mathematical tradition as a whole is much more relaxed and I have seen f(oo), f(-oo) and other shorthands in many places!

@sand15 The Peano curve is not just dense; it passes through every point of the unit square.

You should explain what you are trying to do/know.

eqs:=
int(a*x^2,     x=0..2) = 1/6*8,
int(4 - b*x^2, x=0..c) = 3/6*8,
b*c^2 = 4:
solve([eqs]);

        {a = 1/2, b = 16/9, c = 3/2}

@janhardo Don't forget the constants such as Pi, gamma; they are names, but not variables.

Actually the original mws works after a substitution in the worksheet diag --> Diag
(diag being a protected word after loading the obsolete linalg package). That's all the student did (so, Moise was not quite dead :-)).

@Robert Jantzen 

1. If you (like me!) do not like the new interface (with ribbon tabs) you can still use the old one: simply use maplew-accessibility.exe instead of maplew.exe. I use it in Worksheet mode with 1D input and 2D output.  It is mentained because "Screen Reader Mode"  does not work well in the new interface, see ?accessibility.

2. To see the Executable Group left margin markers  just press the F9 key. I have asked in this forum about the possibilty of reversing this new default setting but I did not receive any answer: https://mapleprimes.com/questions/242029-Hidden-Execution-Group-Boundaries.