Somehow the a~ occurring in the expressions constructed from the use of a label has become a global version. But not always it seems. Here is a simpler version.
You will notice that S and U are identical, but that S2 and U2 are definitely not.
The reason for this behavior with the use of labels I don't know, but clearly it is unfortunate.
I have never used labels myself (well, just now I did).
Also I prefer using assuming instead of assume.

restart;
assume(a>0);
S:=a;
subs(a=2,S);
U:=??; #Insert label (1)
subs(a=2,U);
indets(U) union indets(S);
map(type,%,`local`);
addressof(U);
addressof(S);
about(S);
S2:=a^2;
U2:=??; #Insert label (9)
addressof(S2);
addressof(U2);
indets(U2) union indets(S2);
map(type,%,`local`);
subs(`a~`=2,U2); # The `a~` appearing in the first argument to subs is the global version
sqrt(U2);
sqrt(S2);
about(%);
sqrt(U2) assuming `a~`>0;
about(%);



MaplePrimes13-07-19a.mw

You can use the secant method instead of fsolve.
I also removed option remember and changed the global declaration to local.
The first lines are unaltered, then:
#################
p:=proc(s) local res,F0,F1,F2;
if not type(s,numeric) then return 'procname(_passed)' end if:
res := dsolve({subs(sigma1=s,eq1)=0,eq2=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=0,phi(1-zet)=0,T(0)=0,T(1-zet)=0}, numeric,output=listprocedure);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
evalf((2/(1-zet^2))*Int(F0(eta)*(eta+zet),eta=0..1-zet))-Ree*Pr;
end proc;
plot(p(pp),pp=0..80,adaptive=false,numpoints=25);
s1:=Student:-NumericalAnalysis:-Secant(p(s),s=[60,80],tolerance=1e-8);
res := dsolve({eval(eq1,sigma1=s1)=0,eq2=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=0,phi(1-zet)=0,T(0)=0,T(1-zet)=0}, numeric):
with(plots):
odeplot(res,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1-zet);

In copying your input to MaplePrimes all multiplication signs disappeared. Also pi should be Pi.
I replaced -1.0*Pi by Pi.
I hope I got your equations right.

restart;
eq1:= 9180*Pi - 16295.070*arccos(-1 + 0.001356852103*xu) + 22.1100*sin(-Pi + arccos(-1.0 + 0.001356852103*xu))*(737.0-xu)+7115.0700*arccos(-1.513347022 + 0.002053388090*xu) - 14.6100*sin(-Pi + arccos(-1.513347022 + 0.002053388090*xu))*(737.0-xu) + (0.1050*As*(xu-125)/(1349-xu)) - 0.1050*As = 1688;

eq2:= 9180*Pi- 16295.070*arccos(-1 + 0.001356852103*xu) + 22.1100*sin(-Pi + arccos(-1.0 + 0.001356852103*xu))*(737.0-xu)+7115.0700*arccos(-1.513347022 + 0.002053388090*xu) - 14.6100*sin(-Pi + arccos(-1.513347022 + 0.002053388090*xu))*(737.0-xu)*(0.7370-0.0005*xu) + (0.064260*As*(xu-125)/(1349-xu)) + 0.064260*As= 16014;

plots:-implicitplot([eq1,eq2],As=0..2*10^5,xu=0..1000);
#Using numerical solution:
Sol:=fsolve({eq1,eq2},{xu=400,As=100000});

I didn't get the error you report, but that is surely because my patience is rather limited. I stopped the computation. As far as I can see there is nothing wrong in your input.
I had more luck using a preliminary change of variable and after that fsolve using an initial guess obtained from plotting.

restart;
t0 := arccosh(lambda/(s*(1-lambda)));
f1 := 2*Int(sqrt(lambda-s*cosh(t)/(1+s*cosh(t))), t = 0 .. t0); #The integrand is even!
s := 1/10;
f2:=IntegrationTools:-Change(f1,t=arccosh(x));
plot(f2-Pi/2,lambda=0..1,adaptive=false,numpoints=25);
fsolve(f2-Pi/2=0,lambda=.3);
evalf(eval(f1-Pi/2,lambda=%)); #Check

If you know the values of EL, IN, and q1, you can solve numerically:
res:=dsolve(eval(totsys,{EL=1,IN=1,q1=1}),numeric); #An example
plots:-odeplot(res,[x,w1(x)],0..100);

Replace the direct calls to task with
Start(task,M,1,5);
Start(task,M,2,5);
Start(task,M,5,5);
#See
?Task Programming Model

Note:

Running this several times I get wrong results occasionally. The behavior was the same when I used the full name of Continue (i.e. Threads:-Task:-Continue) inside the procedure 'task'.
I don't know why. I haven't been using Threads before.

From the help page for rsolve:
"The rsolve command can solve linear recurrences with constant coefficients, systems of linear recurrences with constant coefficients, divide and conquer recurrences with constant coefficients, many first order linear recurrences, and some nonlinear first order recurrences."
That rsolve returns unevaluated just means that it cannot solve the problem.
It is easy to see that the problem comes down to solving
eq:= v(n) = v(n-1)*(21*(41/40)^n-v(n-1))/(11*(41/40)^n-v(n-1)) ; #(I converted to rationals)
However,
rsolve({eq,v(0)=1/2},v);
returns unevaluated.
Surely for any finite N you can find v(n) for n <= N, but Maple doesn't find a general formula.
The command
p:=rsolve({eq,v(0)=1/2},v,makeproc);
#or back to floats
#p:=rsolve(evalf({eq,v(0)=1/2}),v,makeproc);
gives you a procedure p s.t. p(n) = v(n), try e.g.
p(4);
evalf(%);

In DEtools there is de2diffop, which was designed for a somewhat more general purpose, but it does give you the characteristic polynomial for a linear homogeneous ode.
Example:
eq:=add(a[i]*diff(x(t),[t$i]),i=0..4)=0;
DEtools[de2diffop](eq,x(t),[lambda,t]);
#The need for the list of two names (here lambda and t), where one name in the present context would appear as enough is due to the broader scope of de2diffop. Try replacing [lambda,t] by [lambda, s], where s is a name. It makes no difference, when the ode has constant coefficients, as should be the case in our present context.

Remark. Clearly you can use de2diffop to make your own simplified version:
charpoly:=proc(eq::{algebraic,equation},yt::function(name),var::name)
   local not_constant_coefficients;
   DEtools[de2diffop](eq,yt,[var,not_constant_coefficients])
end proc;
#The name of the local variable is just for fun, but should alert students to the problem.         
charpoly(eq,x(t),R);
eq2:=add(a[i]*t^i*diff(x(t),[t$i]),i=0..4)=0;
charpoly(eq2,x(t),R);

Are the elements of the outer list either lists of positive integers or lists of sets of positive integers?
In both of your examples that is the case.
Secondly: There are 3 occurrences of "1" in your examples, thus the third "1" should be replaced by "2". But then the set {2,2} is just {2}, but you write {2,3}. Does this mean that you apply the rule again to the result (ignoring that 2 doesn't really occur 3 times because of the {2,2} being {2})?

Anyway, you may have a look at the following idea, which uses StringTools.

restart;
with(StringTools):
L:=[[1],[{1,2}],[{1,2}]];
Lsym:=convert(L,symbol);
S:=indets(L,posint);
St:=map(convert,S,string);
#The following construction allows for checking all the members of St, not only the first (St[1])
try
  n:=SearchAll(St[1],Lsym)[3];
catch:
  n:=0
end try;
#Now only continue if n>0:
Delete(Lsym,n..n);
Insert(%,n-1,convert(S[1]+1,symbol));
parse(%);

 

Did you intend s(b/y) to be s*b/y?
In that case you are missing a multiplication sign and after that correction is made you get 4 solutions in Maple 15 and 17.
If not it is important to know something about the function s.

Maybe this example is helpful:
g[5]:=89;
p:=proc(x) local a; global g;
   a[1]:=14;
   a[2]:=3;
   a[1]*x+a[2]*x^2+g[5]
end proc;
p(y);
p(8);
################
Another example uses the seq modifier:

restart;
g[5]:=89;
p:=proc(X::seq(name)) local a,n,i,x; global g;
   x:=[X];
   n:=nops(x);
   for i to n do a[i]:=i^2 end do;
   add(a[i]*x[i],i=1..n)+g[5]
end proc;
p(y);
p(8);
p();
p(seq(c[i],i=1..6));
?Procedure Parameter Modifiers

There are several issues.

First if you want to change C you need something like
sol:= C1->pdsolve(eval(Eq1,C=C1),ICs union BCs,numeric,time=t,range=0..L,spacestep=.1);
after which you can do e.g.
sol(1):-plot(eta=0,t=0..10);
In the loops you need to change accordingly:
for i from 0 to 600 do
for t1 from 0 to 3 by 1 do
for C1 from 0.1 to 1 by 0.1 do #You cannot start with C1=0: Eq1 is no more a pde of order 2.
p0[C1,t1]:=sol(C1):-value (t=t1):
A[i*0.01]:=p0[C1,t1](i*0.01)
etc.
Notice that A[i*0.01] will end up recording the result corresponding to t1=3 and C1=1 only, which was surely not what you intended since otherwise the two inner loops just waste time!

Remark. You may consider making a change in the time variable letting tau=C*t and v(eta,tau)=u(eta,t):
dchange({t=tau/C,u(eta,t)=v(eta,tau)},Eq1,[tau,v]);
Thereby you have eliminated C from the pde, but it will show up in the boundary condition v(0,tau) = sin(tau/C).

Quoting from the help page for isolate (Maple 12 as it happens):
"Furthermore, whereas isolate returns an equation equivalent to its input, solve returns solutions to its input equations, and can handle systems of equations with multiple solutions."

The last point is relevant here: There are multiple solutions.

eq:=simplify(eq19);
solve(eq,{r});
map(s->s^2,op(%[3]));

After the assignment to g1 I made the following changes:
type(U2y,polynom(anything, y));
#g2:=Int(U2y,x=0..X):
g2:=collect(U2y,y,s->Int(s,x=0..X)): #Int on the coefficients of the polynomial in y.
g22:=(1/delta)*g2:
g3:=U*Ux:
#g33:=Int(g3,x=0..X):
g33:=subs(x=X,U^2/2):
g4:=V*Uy:
type(g4,polynom(anything, y));
#g44:=Int(g4,x=0..X):
g44:=collect(g4,y,s->Int(s,x=0..X)):
P:=subs(x=X,g1)+g22-R*(g33+g44):
type(P,polynom(anything, y));
Eta:=subs(x=X,eta);
P1:=(1/Eta)*int(P,y=0..Eta): #int here, not Int
Po:=subs(X=0,P1):
d1:=(Po-P1):
a0:=subs(k=0.0,d1):
a1:=subs(k=0.1,d1):
a2:=subs(k=0.5,d1):
Then your plotting command gives:

MaplePrimes13-07-10i.mw

You should not use linalg. It was replaced by LinearAlgebra in Maple 6 (yes that is a long time ago).
Is this what you want:
A := <`&varphi;`[1, 1]|`&varphi;`[1, 2]|`&varphi;`[1, 3]|`&varphi;`[1, 4]>;
Notice the vertical bars.
Another syntax is:
A:=Matrix([`&varphi;`[1, 1],`&varphi;`[1, 2],`&varphi;`[1, 3],`&varphi;`[1, 4]]);