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;

restart;
sys := {6*(diff(a(t), t))^2+12*a(t)*(diff(a(t), t$2))-3*a(t)^2*phi(t)^(-2*c)*sqrt(1-alpha*(diff(phi(t), t))^2), 2*c*a(t)^3*phi(t)^(-2*c-1)*sqrt(1-alpha*(diff(phi(t), t))^2)-3*alpha*a(t)^2*phi(t)^(-2*c)*(diff(a(t), t))*(diff(phi(t), t))/sqrt(1-alpha*(diff(phi(t), t))^2)-alpha*a(t)^3*phi(t)^(-2*c)*(diff(phi(t), t$2))/sqrt(1-alpha*(diff(phi(t), t))^2)+2*c*alpha*a(t)^3*phi(t)^(-2*c-1)*(diff(phi(t), t))^2/sqrt(1-alpha*(diff(phi(t), t))^2)-alpha^2*a(t)^3*phi(t)^(-2*c)*(diff(phi(t), t))^2*(diff(phi(t), t$2))/(1-alpha*(diff(phi(t), t))^2)^(3/2)};
R(t) := 6*((diff(a(t), t))^2/a(t)^2+(diff(a(t), t$2))/a(t));
W(t) := -phi(t)^(-2*c)*sqrt(1-alpha*(diff(phi(t), t))^2)/(1/a(t)^3+a(t)^3+phi(t)^(-2*c)/sqrt(1-alpha*(diff(phi(t), t))^2));
inc:=a(0)=1.5,D(a)(0)=0.1,phi(0)=1,D(phi)(0)=-0.5;
sol:=dsolve({sys[], inc}, numeric, parameters=[c, alpha]);

sol(parameters=[1,1]);
plots:-odeplot(sol, [[t,a(t)], [t,R(t)]], t=0..10, color=[red,blue]);
 

I guess that  :=  (assignment sign) should be instead of  =  :

restart; 
M := 5;
for i while i <= M do N[i, 0](u) := 1; N[0, i](u) := 0 end do;

Another way is the usage  assign  command for a multiple assignment (without any loops):

 assign(seq(op([N[i, 0](u)=1, N[0, i](u)=0]), i=1..M));
 

restart;
f:=(x,c,k)->c*x/(k+x):
plot([f(x,5,2), f(x,5,3)], x=0..50, color=[red,blue], legend=["c=5 \n k=2","c=5 \n k=3"], size=[800,400]);

                       

The corrected version:

1_(1)_new.mw

restart;
eqns := {(1-sin(b_u)*sin(s)/(cos(b_u-y_f)*cos(s-y_f)))*(1/cos(s-y_f)^2-1/sin(s)^2) = -(cot(s)+tan(s-y_f)+Z)*sin(b_u)*(cos(s)/cos(s-y_f)+sin(s)*sin(s-y_f)/cos(s-y_f)^2)/cos(b_u-y_f)};
b_u := 1/tan(0.8);
Z := 892/(27417000*f_z);
y_f := 9*Pi/180;
S:=[seq(rhs(fsolve(eval(eqns, f_z=a), s)[]), a=0.00005..0.0005, 0.000001)];

S is the list of desired solutions.

Addition. The plot:

plot([seq(a, a=0.00005..0.0005, 0.000001)], S);