You may obtain animations using Explore.
Here is a simple example, similar to your problem.

restart;
A:=proc(c)
local p1,p2,d,x,t ;
p1:=plot(1 + c*sin(t), t=0..2*Pi, coords=polar);
d:=fsolve(x^5+x = c^5+1, x);
p2:=plot(2*x+d, x=-2..1);
plots:-display(p1,p2, view=[-2.5 .. 2.5, -1..4]) 
end proc:

Explore(A(c), c=0..2.6, animate);

BISSEC := proc (PP, UU, VV) 
uses LinearAlgebra,plots;
local P:=PP,U:=UU, V:=VV,
a, b, eq1, eq2, M1, M2, t, PU, PV, bissec1, bissec2; 
a := (P-U)/Norm(P-U, 2)+(P-V)/Norm(P-V, 2); M1 := P+a*t; b := (P-U)/Norm(P-U, 2)-(P-V)/Norm(P-V, 2); M2 := P+b*t; 
eq1 := op(eliminate({x = M1[1], y = M1[2]}, t)); 
eq2 := op(eliminate({x = M2[1], y = M2[2]}, t)); 
P := convert(P, list); U := convert(U, list); 
V := convert(V, list);
PU := plot([P, U]); PV := plot([P, V]); bissec1 := implicitplot(op(eq1[2]), x = 0 .. 5, y = 0 .. 10, color = red); bissec2 := implicitplot(op(eq2[2]), x = 0 .. 5, y = 0 .. 10, color = green); display([bissec1, bissec2, PU, PV], scaling = constrained) end proc:

A := `<,>`(4, 8): B := `<,>`(4, 2):; C := `<,>`(1, 4):
#debug(BISSEC);
BISSEC(A, B, C);

P:=n -> binomial(n-1,3) - binomial(floor(n/2),3);

Usually, select or map produce shorter and faster solutions than a loop.

select(`>`, L, 0) /~ 2;

DD:=(n,alpha) -> Matrix(n+1, (i,k) -> sum(omega(i-1,j,k-1), j=ceil(alpha)..n)):
DD(2,alpha);

It remains to define your function omega:

omega := (i,j,l) -> ...

Yes, in Maple 2019 it is recommended to use

g := x-> local j; sum(..., j=1..5);

I have used DirectSearch, a reliable and free global optimizer (see  Maple Application Center).
Note that your model does not seem to be perfect for your data.

plots:-display(plot(F, t=0..500), plot(X,Y, style=point, color=blue, symbolsize=16));

The 0 value can be obtained as an Abel (generalized) summation:

sum(sin(Pi*n/2)*sin(n*Pi*(x + 1)/2)*cos(n*Pi*t/2) * z^n, n=1..infinity) assuming z>0,z<1:

limit(%, z=1, left);
      0

Your piecewise function f is not a "classical" solution of the ODE y' = 0  because as you see, f is not differentiable at some points.

There are also "generalized" solutions,  see wiki, where the function must satisfy the ODE only a.e. (almost everywhere, i.e. except on a null set), but the function must be absolutely continuous. Your piecewise funcion is not even continuous, so it is still not a solution in this wider sense.
Note that there are other contexts (e.g. distributional derivatives), but I'll stop here.

You essentially compare

eval(floor(n), n=1);   # 1
limit(floor(n), n=1);  # undefined

and conclude 1 = undefined (+ something about a prize).

You should be aware that floor is not continuous at 1.

LibLocation:="C:/temp/MyLib.mla": # dir with write access
LibraryTools:-Create(LibLocation); # if the library is new
LibraryTools:-Save(MyMat, LibLocation);
restart;  # check
LibLocation:="C:/temp/MyLib.mla";
libname:=LibLocation, libname;
with(MyMat);

(The path is for Windows.)

The command works for any m, provided that A in invertible mod m. You can even use e.g.

A:=<2,7;3,2>:
1/A mod 10;
   

% . A mod 10;  # check

You may also use:
LinearAlgebra:-Modular:-Inverse(m, A);
[which is faster for large matrices]

Note that A is invertible mod m iff det(A) and m  are relatively prime.

Maple is not needed for this. k=0 obviously satisfies the condition.
For any integer k>0, a(n) cannot be integer if n is large enough (because the distance between two consecutive squares --> oo).

Edit. The answer refers to the recurrence  a(n+1) = 3*a(n)+sqrt(4*a(n)^2+k), a(1)=1.
As I see, the OP changes (without notice) the problem every  half hour.
Answering the new one is very risky!

F := unapply(int(z^(-1/2), z), z);
              F := z -> 2*sqrt(z)
limit(F(-1-I*y),y=0, right) - F(1);
             -2 - 2*I

If you know the factorization, just add it to the remember table of ifactors:

ifactors(n):=[1, [[p, 1], [q, 1]]];

and then call ModularSquareRoot.