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);

 

1. The 1000-th prime is  7919  not  1979

2. Use  unapply  command instead of  -> :

restart;
P:= 7919; # (this is the1000-th prime)
S:=5000; k:=5:
a:=[0$i=1..k];
f:=[0$i=1..k];
a[1]:=S;
f:=unapply(a[1], x);

with(RandomTools):
for i from 2 to k do
a[i]:=Generate(integer(range=1..P-1));
f:=unapply(f(x)+a[i]*x^(k-1), x);
print(f(x));
od;