@Carl Love  Here is it:

Download permanouuuuuuuuuuuuuuuuuun_new.mw

@JAMET  To get the curve equation from equations  t  and  t1  in the form  F(x,y)=0  you just need to eliminate the parameter  m   from these equations:

restart;  
t := y-m*x-p/(2*m) = 0; 
t1 := y+x/m+(1/2)*p1*m = 0; 
Sol:=eliminate({t,t1}, m);
Eq:=op(Sol[2])=0;
plots:-implicitplot(eval(Eq,[p=1,p1=-2]), x=-2..2,y=-2..2,color=red, scaling=constrained, size=[200,600]); 

@vv  I don’t understand where the discontinuous points come from? The integrand is obviously continuous. Isn't its antiderivative function to be continuous?

For example, if we plot  F , we see a continuous curve:

F:=Int(sqrt(1+((k*Pi)/l*cos((Pi*x)/l))^2), x=0..X):
plot(eval(F,[k=2,l=2]), X=0..6);

@Mariusz Iwaniuk  Of course you're right. The integrand  sqrt(1 + k^2*Pi^2*cos(Pi*x/l)^2/l^2)  in the original integral is positive, therefore the definite integral  (if  b>a) is also positive. Maple seems to find this primitive function  F1 =  -csgn(sin(Pi*x/l))*Elliptic E(cos(Pi*x/l), I*k*Pi/l)*l/Pi  incorrectly.

@Carl Love  for this info. I fixed the procedure  P  code with this in mind.

@ActiveUser  Present here in text form code that does not work in Maple 12 (so that I can test it). The upper value for the variable  n  in the procedure  P  is not indicated , because it is unknown in advance and can be large enough if the numbers  a  and  b  are close.

PS. I slightly modified the procedure  P  code. Check again how it works.

@Stretto   I did not study all these options in detail, but I just offered you another way in which you first find all these points, and then consciously delete them. Of course, you yourself must choose the approach that you like best.

As for me, I always try to choose an approach where I clearly understand what is happening and it is easy for me to manage what is happening.

@Carl Love  Your code does not work in Maple 2018.2 and in older versions, at least since Maple 2015.

@Carl Love  Have you checked if this works?

restart;
r:= t->[t, t^2, t^3];
a:=1;
D(r)(a)*~(t-a)+~r(a);

@ActiveUser  If I understand correctly, then you can just take the arithmetic mean of any neighboring numbers.
For your example:
 

restart; 
sourcesamples :=[evalf(-1-sqrt(7)), -2, 1, evalf(-1 + sqrt(7)), 2];
n:=nops(sourcesamples);
samples:=sort(convert({ceil(sourcesamples[1])-1, seq(op([sourcesamples[i],(sourcesamples[i]+sourcesamples[i+1])/2,sourcesamples[i+1]]), i=1..n-1),floor(sourcesamples[n])+1},list),(x,y)->evalf(x)<evalf(y));

@ActiveUser  I am not familiar with all this issue.

@ActiveUser 

1. You can use the  op  command for this: 

Sol:=solve(a*x^2 + b*x + c=0, x, parametric=true, real);
op(Sol);
op(9,Sol);
op(9..10,Sol);

2. I do not know the commands from this package and have never used it.

@Carl Love I do not understand what your doubts are based on. As for the work of the value procedure, I would call it not an error, but simply an unsuccessful design. The essence of the problem is that this procedure does not calculate any derivatives (but only the values of the functions themselves from the system), and returns  0  for any derivatives. Of course, it would be better if in the case of derivatives, it returns the input unevaluated or returns an error message. I believe that spacestep and timestep options are not connected in any way with this problem.

Compare the angular coefficients of the tangents at  y = 0  with the numerical results (this to some extent confirms the correctness of the results):

restart;
a := 0.7: L := 8: HAA := [0, 2, 5, 10]: h:=0.001:
PDE1 := diff(u(y, t), t) = diff(u(y, t), y, y)+diff(diff(u(y, t), y, y), t)-b*u(y, t)+T(y, t):
PDE2 := diff(T(y, t), t) = (1+(1+(a-1)*T(y, t))^3)*(diff(T(y, t), y, y))+(a-1)*(1+(a-1)*T(y, t))^2*(diff(T(y, t), y))^2+T(y, t)*(diff(T(y, t), y, y))+(diff(T(y, t), y))^2:
ICandBC := {T(L, t) = 0, T(y, 0) = 0, u(0, t) = t, u(L, t) = 0, u(y, 0) = 0, (D[1](T))(0, t) = -1}:
for i from 1 to nops(HAA) do
b := HAA[i];
pds[i] := pdsolve({PDE1,PDE2}, ICandBC, numeric);
f:=(y0,t0)->rhs((pds[i]:-value(u(y,t),t=t0)(y0))[3]);
A[i]:=(f(h,1)-f(0,1))/h;
end do:
seq(A[i], i=1..4);
plots:-display(< `<|>`(seq(pds[i]:-plot(u(y,t), y=0..3, t=1, scaling=constrained),i=1..4))>);

If you need a procedure, then do so (I just rewrote vv's code as a procedure):

Squares:=(M,N)->plots:-display(seq(seq(plottools:-polygon([[m,n],[m+1/2,n+1/2],[m+1,n],[m+1/2,n-1/2]], color=`if`(m::odd,red,blue)), m=0..M),n=0..N), axes=none, scaling=constrained):

Examples of use:
Squares(2,2);
Squares(5,5);
Squares(3,2);

@vv 

sys := [x^2 + y^2 - x*y - 1 = 0, y^2 + z^2 - y*z - a^2 = 0, z^2 + x^2 - x*z - b^2 = 0]:
ab:={a=10,b=104/10}:
solve(eval(sys,ab), explicit):
evalf(%);

           {x = 6.676479229-2.966007548*I, y = 5.939434090+4.226624806*I, z = 12.50085790+.1379423560*I}, {x = -6.676479229+2.966007548*I, y = -5.939434090-4.226624806*I, z = -12.50085790-.1379423560*I}, {x = 6.676479229+2.966007548*I, y = 5.939434090-4.226624806*I, z = 12.50085790-.1379423560*I}, {x = -6.676479229-2.966007548*I, y = -5.939434090+4.226624806*I, z = -12.50085790+.1379423560*I}, {x = .1989718373, y = 1.084528494, z = 10.49807685}, {x = -.1989718373, y = -1.084528494, z = -10.49807685}, {x = 1.141952419, y = .4227976407, z = -9.781889234}, {x = -1.141952419, y = -.4227976407, z = 9.781889234}

Or even simpler:

sys := [x^2 + y^2 - x*y - 1 = 0, y^2 + z^2 - y*z - a^2 = 0, z^2 + x^2 - x*z - b^2 = 0]:
ab:={a=10,b=10.4}:
solve(eval(sys,ab));

Output is the same.