From your question, I first learned about  Interpolation  package. But when I tried to deal with his work when working with your examples, the message always appears "Kernel connection has been lost"

All your examples are easy to solve in  CurveFitting  package. Only it is necessary to use an exact arithmetic, otherwise due to rounding errors Maple writes that at the node points  (x=0,1, .. ,10))  the derivative is not defined.

restart;
with(CurveFitting):
points := [seq(x, x = 0 .. 10)]; 
data := [seq(sqrt(1+x^2), x = 0 .. 10)];
f:=unapply(Spline(points, data, x), x);
int(f(x), x = 0 .. 1);
evalf(%); 
eval(diff(f(x), x), x = 2);
evalf(%);
plot(f(x), x = 0 .. 10);
plot([f(x), diff(f(x), x)], x = 0 .. 10);
 

Because  I  is the imaginary unit in Maple, then I use  T  instead. Below we build a table  T  that allows us to find the values  T[i1,i2,i3,i4]  for the indices  [i1,i2,i3,i4] . I think that this is something equivalent to a tensor.

restart;
for i1 from 1 to 3 do
  for i2 from 1 to 3 do
    for i3 from 1 to 3 do
      for i4 from 1 to 3 do
        if i1=i2 and i2=i3 and i3=i4 then
        T[i1,i2,i3,i4]:=1.0;
        else
        T[i1,i2,i3,i4]:=0.0 end if;
        end do:
      end do:
    end do:
  end do:
T:=convert(T, table);
                

Examples of use:

T[1,2,3,4];
T[2,2,2,2];
                                           0.
                                          1.

Of course, we can define this much shorter as a function of 4 variables:

T:=(i1,i2,i3,i4)->`if`(i1=i2 and i2=i3 and i3=i4, 1., 0.):

Examples of use:

T(1,2,3,4);
T(2,2,2,2);

I replaced  sum  with  add  and got the following result:

restart;
M:=4: N:=2: alpha:=1:
 
add(add(((-1)^i2*GAMMA(N-i2+alpha)*2^(N-2*i2)/(GAMMA(alpha)*factorial(i2)*factorial(N-2*i2)*(N-2*i2+1))*(GAMMA(k+1)*(k+alpha)*GAMMA(alpha)^2/(Pi*2^(1-2*alpha)*GAMMA(k+2*alpha))))*(add((1/2)*(-1)^i*GAMMA(k-i+alpha)*2^(k-2*i)*(1+(-1)^(N-2*i2+1+k-2*i))*GAMMA((1/2)*N-i2+1+(1/2)*k-i)*GAMMA(alpha+1/2)*L[k]/(GAMMA(alpha)*factorial(i)*factorial(k-2*i)*GAMMA(alpha+3/2+(1/2)*N-i2+(1/2)*k-i)), i = 0 .. floor((1/2)*k))), i2 = 0 .. floor((1/2)*N)), k = 0 .. M);

                                               -(1/6)*L[1]+(1/6)*L[3]
                                              

You can easily calculate this using  LinearAlgebra  package (here we have a more compact syntax than the one suggested above):

P1 := [-17/12,11/36,-65/36]:
P5 := [-11/12,-7/36,-47/36]:
P6 := [-2/3,-1/9,-8/9]:
P8 := [-5/4,-1/12,-17/12]:
V1:=convert(P8-P5, Vector);
V2:=convert(P6-P5, Vector);
V3:=convert(P1-P5, Vector);
with(LinearAlgebra):
V1 &x V2 . V3;
 

restart;
A:=proc(x::list)
local N, L;
N:=nops(x);
L:={$1..N};
Matrix(N,(i,j)->`if`(i<>j,mul(x[i]-x[k],k=L minus {i,j})/mul(x[j]-x[k],k=L minus {j}), add(1/(x[i]-x[k]), k=L minus {i})));
end proc:

Example of use:

A([a,b,c,d]);
           

I've adjusted your syntax a little. Maple finds the exact solution of this equation in a closed form:

restart;
sys[1] := [-(diff(u(x, t), t, t))-(diff(u(x, t), x, x))+u(x, t) = 2*exp(-t)*(x-(1/2)*x^2+(1/2)*t-1), u(x, 0) = x^2-2*x, u(x, 1) = u(x, 1/2)+((1/2)*x^2-x)*exp(-1)-((3/4)*x^2-(3/2)*x)*exp(-1/2), u(0, t) = 0, D[1](u)(1, t) = 0]:
ans:=pdsolve(sys[1]);
pdetest(ans, sys[1]);
plot3d(eval(u(x,t), ans), x=0..2, t=0..2);

CalculesS:=proc(n)
local S1, S2, i, j;
S2:=0:
for i from 1 to n do
S1:=1;
for j from i to n do
S1:=S1*1/j^i;
od;
S2:=S2+S1;
od:
S2;
end proc:

Example of use:

CalculesS(100);

Of course, it's easier to write without any loops, if you use the commands  product  and  sum :

S:=n->sum(product(1/j^i, j=i..n), i=1..n); 
S(100);
 

I did not find any contradictions. Of course, ans1  and  ans2  differ in form, but in fact they determine the same set of points (the intersection of two circles). ans1  and  ans3  are absolutely identical:

is(convert~(convert(ans3, set), set)={ans1});
                                                                  true
 

Maple calculates the values of this function in the range  -1..1 as complex. Look

seq(LegendreQ((1/2)*sqrt(5)-1/2, x), x=-0.9..0.9, 0.1);

                -1.149108610+1.159683101*I, -1.227886365+.7316329045*I, -1.247731828+.4486859813*I, -1.240765970+.2275713122*I, -1.217080008+0.4156209706e-1*I, -1.181045854-.1215055276*I, -1.134806859-.2682523506*I, -1.079402624-.4027068739*I, -1.015207630-.5275156677*I, -.9421078357-.6445179356*I, -.8595474579-.7550479980*I, -.7664842691-.8601073870*I, -.6612453970-.9604688726*I, -.5412284504-1.056742479*I, -.4023020093-1.149419091*I, -.2375185873-1.238900244*I, -0.3395679811e-1-1.325519053*I, .2369504060-1.409555298*I, .6644541489-1.491246518*I

You can use  InertForm  package for this:

InertForm:-Parse("x+(y^(2+x)-4)/3");
lprint(%);
value(%);

Curves:=plot([tan(x), tan(Pi/4)+D(tan)(Pi/4)*(x-Pi/4), arctan(x), arctan(1)+D(arctan)(1)*(x-1)], x=0..3, 0..3, color=[red,yellow,blue,green], discont, legend=[tan(x),"Tangent for tan(x)", arctan(x),"Tangent for arctan(x)"]):
Points:=plot([[[Pi/4,1]], [[1,Pi/4]]], style=point, color=[red,blue], symbol=solidcircle, symbolsize=10):
plots:-display(Curves, Points, scaling=constrained, size=[400,400]);

BVP := [4*(diff(u(x, t), t))-9*(diff(u(x, t), x, x))-5*u(x, t) = 0, u(0, t) = 0, u(6, t) = 0, u(x, 0) = sin((1/6)*Pi*x)^2]:
sol:=pdsolve(BVP):        
U := unapply(rhs(sol), [x, t]);  
plot3d(add(op(1,U(x, t)), n=1..20, 2), x = 0 .. 6, t = 0 .. 4);

To find the derivative of the function  f  with respect to the first variable for certain values of these variables, you must first define the function  f  itself in terms of its variables, for example:

restart;
f:=(u,v)->u^2+u*v;
D[1](f)(z*sqrt(a/nu[f]), U*t/(2*x));

# The first way
seq(seq(exp(a*Pi/b), b=1..10), a=1..15);
nops([%]); 

# The second way
k:=0:
for a from 1 to 15 do
for b from 1 to 10 do
k:=k+1; L[k]:=exp(a*Pi/b);
od: od:
convert(L, list)[];  

The results are identical.

Place:=proc(L,n)
local k;
uses ListTools;
k:=Search(n,L);
if k<>0 then return k else NULL fi;
end proc:

Examples of use:

Place([1,2,4,5], 4);
Place([1,2,4,5], 6);