DS:=proc(X::set, C::set(set)) # Dynkin system
local D:=C union {X}, S, T, ok:=false;
if not `union`(C[]) subset X then error "C must contain subsets of X" fi;
while not ok do
  ok:=true;
  for S in D do for T in D do
    if S subset T and not (T minus S) in D then D:= D union {T minus S}; ok:=false fi;
  od od
od;
D
end proc:

DS({1,2,3}, {{1},{2}});nops(%);
     {{}, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}

                                8

DS({1,2,3,4,5}, {{1,4},{2,4}, {1}}); nops(%);
{{}, {1}, {2}, {4}, {1, 2}, {1, 4}, {2, 4}, {3, 5}, {1, 2, 4},  {1, 3, 5}, {2, 3, 5}, {3, 4, 5}, {1, 2, 3, 5}, {1, 3, 4, 5},  {2, 3, 4, 5}, {1, 2, 3, 4, 5}}

                               16

DS({1,2,3,4,5}, {{2,4,5},{1,4,5}}); nops(%);
     {{}, {1, 3}, {2, 3}, {1, 4, 5}, {2, 4, 5}, {1, 2, 3, 4, 5}}

                               6

Download geometry-parab.mw

You probably want non-collinear vectors.

restart;
P:=(p,q,r) -> [p^2-q^2-r^2 ,2*p*q, 2*p*r]:
n:=13:
pyt:=[seq](seq(seq(P(p,q,r), p=floor(sqrt(q^2+r^2))+1..n), q=1..n),r=1..n):
pyt:=select( u -> evalb(igcd(u[])=1), pyt):
for i to nops(pyt) do for j from i+1 to nops(pyt) do
 A:=pyt[i]; B:=pyt[j];
 N2:=(A[2]*B[3]-A[3]*B[2])^2+(A[1]*B[3]-A[3]*B[1])^2+(A[1]*B[2]-A[2]*B[1])^2;
 if issqr(N2) then lprint(A,B) fi
od od:

[31, 24, 12], [3, 24, 16]
[19, 48, 12], [3, 16, 24]
[19, 48, 12], [3, 96, 80]
[27, 96, 16], [13, 144, 36]
[63, 120, 20], [115, 120, 48]
[63, 120, 20], [7, 24, 60]
[107, 144, 24], [3, 16, 24]
[31, 12, 24], [3, 16, 24]
[3, 24, 16], [19, 12, 48]
[3, 24, 16], [107, 24, 144]
[131, 72, 48], [7, 24, 60]
[61, 72, 36], [83, 120, 144]
[7, 60, 24], [9, 84, 28]
[7, 60, 24], [131, 48, 72]
[7, 60, 24], [103, 96, 120]
[7, 60, 24], [63, 20, 120]
[47, 140, 40], [5, 260, 208]
[13, 144, 36], [3, 80, 96]
[11, 48, 36], [39, 64, 48]
[11, 48, 36], [119, 96, 72]
[11, 48, 36], [21, 176, 132]
[11, 48, 36], [69, 208, 156]
[39, 64, 48], [119, 96, 72]
[39, 64, 48], [21, 176, 132]
[39, 64, 48], [69, 208, 156]
[119, 96, 72], [21, 176, 132]
[119, 96, 72], [69, 208, 156]
[27, 160, 60], [3, 96, 80]
[19, 12, 48], [3, 80, 96]
[61, 36, 72], [83, 144, 120]
[11, 36, 48], [39, 48, 64]
[11, 36, 48], [119, 72, 96]
[11, 36, 48], [21, 132, 176]
[11, 36, 48], [69, 156, 208]
[39, 48, 64], [119, 72, 96]
[39, 48, 64], [21, 132, 176]
[39, 48, 64], [69, 156, 208]
[119, 72, 96], [21, 132, 176]
[119, 72, 96], [69, 156, 208]
[49, 72, 72], [83, 144, 120]
[49, 72, 72], [83, 120, 144]
[103, 120, 96], [7, 24, 60]
[103, 120, 96], [3, 80, 96]
[69, 132, 88], [33, 260, 156]
[7, 24, 60], [9, 28, 84]
[115, 48, 120], [63, 20, 120]
[103, 96, 120], [3, 96, 80]
[3, 96, 80], [13, 36, 144]
[27, 16, 96], [13, 36, 144]
[69, 88, 132], [33, 156, 260]
[3, 80, 96], [27, 60, 160]
[21, 176, 132], [69, 208, 156]
[47, 40, 140], [5, 208, 260]
[21, 132, 176], [69, 156, 208]
 

Download ALG.mw

@Zeineb Maple can help to compute it, but only if one uses the right theorem.
 

restart;
alpha := x->floor(x):
f:=x->x^2+1:
int(f(x)*diff(alpha(x),x), x=0..3);  # wrong
#                               0
f(3)*alpha(3)-f(0)*alpha(0) - int(alpha(x)*diff(f(x),x), x=0..3); # correct
#                               17

As I have mentioned, Int_[0..3] alpha d alpha   does not exisist.

If P[1], ..., P[4] are your points, consider a fifth point P[5] of unknown coordinates, say (a,b).
Then take the conic (via its quadratic equation) thru P[1], ..., P[5] and impose the condition to be a parabola.
In general there will two parabolas, but it is also possible to have 0 (!). E.g. a parabola cannot pass thru the vertices of a square.

Student:-Basics:-LongDivision(A,B,x, output=printf);

                4/3*x^2+80/3*x+1552/3
              -----------------------
3*x^2-60*x+36 ) 4*x^4
                4*x^4-80*x^3+48*x^2
                -------------------
                     80*x^3-48*x^2
                     80*x^3-1600*x^2+960*x
                     ---------------------
                           1552*x^2-960*x
                           1552*x^2-31040*x+18624
                           ----------------------
                                 30080*x-18624
 

You must export in a directory with write permission. By default the current directory is Maple's directory (without such permission).
So, you must change the directory, e.g. use "c:/temp/Test.jpg", (if you have such a directory),
or, execute first
currentdir(kernelopts(homedir));

Your second equation Y, equivalent to Y = 0,  reduces to a = 0, and does not depend on m. So, you have nothing to eliminate.

For higher-order functionals, use variables to represent derivatives. E.g.
L1:=eval(L, [diff(x(t),t,t)=diff(x1(t),t), diff(theta(t),t,t)=diff(theta1(t),t)]):
L2:=L1 + A*(diff(x(t),t) - x1(t))^2 + B*(diff(theta(t),t) - theta1(t))^2:
EulerLagrange(L2, t, [x(t), theta(t), x1(t), theta1(t)]);
 

Such an assignment is not possible due to the standard evaluation in any CAS: psi(t) would be evalaluated to psi(t)(t)(t)... (i.e. a nonsense).
It is possible to use macro (for input only)  or alias (for both input and output). You must be aware that in both cases psi(a) will be evaluated to psi(t)(a).

# x,y,z = the unknowns, top to bottom
# h1,h2,h3  = halves of the heights 
# a,b,c,d,e = the bases of the triangles
restart
eq:=
[260=a*h2, 247=b*h2, 287=d*h3, 209=e*h3,
x=a*h1, y=c*h2, z=(b+c)*h3,
h1/(h1+h2)=a/(b+c), h1/(h1+h2+h3)=a/(d+e)]:
S:=eliminate(eq,{h1,h2,h3,a,b,c,d,e}):
elim:=S[1]:
XY:=S[-1]:
XY1:=[isolve(XY[1])]:
xy:=select( u -> (eval(x,u)::posint and eval(y,u)::posint), XY1):
for v in xy do
  zz:={solve( eval(XY[2],v) )};
  for z1 in zz do  if z1::posint then print(v[],'z'=z1, eval(elim,[v[],'z'=z1])) fi od;
od:

So, two solutions found.
I don't see any relation with the star (100). Maybe I miss something.

It often works.

restart;
r := x*(diff(theta(t), t))^2+y*(diff(varphi(t), t))^2+z*(diff(theta(t), t$2))+w*diff(varphi(t), t$2):
g := (4*(f+T))*(diff(theta(t), t))^2+u*(diff(varphi(t), t))^2+(f+9)*(diff(theta(t), t$2))+4*s*diff(varphi(t), t$2):

eval(r-g, [theta=(t->t^3), varphi=(t->t^7)]): # or similar
solve([coeffs(%, t)], [x,y,z,w]);
simplify(eval(r-g, %[1]));   #check
#          [[x = 4*f+4*T, y = u, z = f+9, w = 4*s]]
#                               0