To solve this problem you can use  PosIntSolve  procedure from this post. This procedure allows you to find the number of non-negative (or positive) integer solutions of a linear Diophantine equation with positive coefficients. For convenience, here are the text of the procedure and 2 examples.

restart;

PosIntSolve:=proc(L::list(posint), n::nonnegint, s::symbol:=nonneg)
local N, n0;
option remember;
if s=pos then n0:=n-`+`(op(L)) else n0:=n fi;
N:=nops(L);
if N=1 then if irem(n0,L[1])=0 then return 1 else return 0 fi; fi;
add(PosIntSolve(subsop(1=NULL,L),n0-k*L[1]), k=0..floor(n0/L[1]));
end proc:

The problem: find the number of non-negative integer solutions of equations  51*x1+53*x2+50*x3=197  and   51*x1+53*x2+50*x3=2697  .

Solutions:

PosIntSolve([51, 53, 50], 197);
PosIntSolve([51, 53, 50], 2697);
                                                             0
                                                             28
 

Replace  "5"  temporarily some symbol:

applyrule(a*x+a=a*``(x+1), subs(5=a, sqrt(5*x+5+y))):
subs(a=5, %);
                             sqrt(5*(x+1)+y)

Actually no error. When you differentiate the internal function  r*sin(theta)*(r*cos(theta)-a)^(-1)  by a-variable, then the factor (-1) occurs 2 times.

The result above is incorrect, since the straight line and the parabola in this example intersect not at one, but at two points. To avoid errors for different curves, use different parameter names. Here is the right solution:

Download parametrics.mw

Red arrows (tangent vectors) show the direction of increasing parameter:
 

Download arrows1.mw

Try the following:

plot([x+deltax, z, x=0..l]);
 

Here is the plot that OP wished according to his link:

restart;
f:=x->sqrt(1-x^2):
g:=x->-sqrt(1-x^2):
y1:=f(a)+D(f)(a)*(x-a):
y2:=g(a)+D(g)(a)*(x-a):
y3:=[-1,t,t=-1.5..1.5]:
y4:=[1,t,t=-1.5..1.5]:
Opt:=x=-2..2,color="Blue":
plots:-display(plot([seq(op([y1,y2]),a=-0.99..0.99,0.02)],Opt),plot([y3,y4],Opt), view=[-2..2,-1.5..1.5], scaling=constrained, size=[700,500],tickmarks=[spacing(0.5),spacing(0.5)],  axis[1]=[gridlines=40], axis[2]=[gridlines=30]);

Output:
                    

L := [8,7,9,12];
L1:=[seq([i, L[i]], i=1..nops(L))];
map(t->t[1], sort(L1, (x,y)->y[2]<x[2]));
 

Your second result obtained with erroneous syntax (numeric option is  misplaced) is incorrect. To test this, I did both calculations for  T=1, after making some simplifications to your commands (Pi should be instead of pi, and also removed some commands and options that do not affect the result). I don't know how Maple could get it but the second result is also erroneous because the integral in which the integrand is positive must also be positive.

Download int_(2)_new.mw

What function do you call monotonically increasing (or monotonically decreasing) on some multidimensional set? The issue is that there is a natural order on the number axis: for any two different points, we know for sure which one is greater. And in multidimensional space? For example, on 2D which point is greater: (0,1) or (1,0)?

I do not understand why such a complicated code (nested loops and etc.) for a simple task. You can immediately build all the necessary plots, if you know their equations, and then get the data using  plottools:-getdata  command:

restart;  
H:= 3+2*tanh(x-6);
P:=plot([H, -H, 0], x=0..10, linestyle=[1,1,3], color=black, thickness=[2,2,0], view=[-1..12,-6..6], axes=box);
interface(rtablesize=1000);
map(t->t[3],[plottools:-getdata(P)]);

restart;
f:=ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2));
applyop(t->numer(t)/expand(denom(t)), 1, f);
simplify(%);
applyrule(ln((-1+x)^2)=2*ln(abs(x-1)), %);
sort(%);
                

Edit.
                    

I do not know of any way to solve these examples in Maple in the above formulation. The first 2 examples can be easily solved using  isolve  command to solve equations in integers:
 

We see that the first assertion is true for any integer greater than or equal to 3, and the second one is true for an infinite set of integers depending on the integer parameter  _Z1

Download prop.mw

Perhaps you mean Puiseux series, which, for example, allow locally to obtain individual branches of algebraic functions near specific points. These expansions may have fractional exponents.

Example:

with(algcurves):
f:=x^3+y^3-3*x*y;
plots:-implicitplot(f, x = -3 .. 3, y = -3 .. 3, gridrefine = 5);
puiseux(f, x = 0, y, 6);
                          

See the quality of the approximation on the plot:
f:=x^3+y^3-3*x*y;
A:=plots:-implicitplot(f, x = -3 .. 3, y = -3 .. 3, gridrefine = 5):
P:=algcurves:-puiseux(f, x = 0, y, 6);
B:=plot([P[]], x = -2 .. 2,color=[red,blue], thickness=2):
plots:-display(A,B, scaling=constrained);
 

Also see help on  ?examples, algcurve

R:=z-> 1-cos(Pi*z)^2;
for z from 0 to 1 by 0.1 do
R(z);
od;