@digerdiga Only solutions analytic at 0 are found, but you can do like this when r = 3/2:
eval(ode,X(x)=x^(3/2)*u(x));
ode1:=expand(%/x^(3/2));
diffeqtorec(ode1, u(x), a(n));

for n from 1 to 8 do
  if n=3 then next end if;
  print(n)
end do;

for n from 1 to 8 do
  if n=3 then next end if;
  print(n)
end do;

@snijman The solution to the system eq1,eq2 is the intersection point of the two curves made by implicitplot. So we pretty much know the answer from the plot. This point is then used as an initial point to feed to fsolve. And yes, fsolve may not succeed if not provided with either ranges or initial points, or it may take much longer, or it may find another solution than the one intended.
If you have no idea in which domain to search for a solution the first could be to use fsolve without initial point or ranges. Otherwise try plotting in various ways.

@snijman The solution to the system eq1,eq2 is the intersection point of the two curves made by implicitplot. So we pretty much know the answer from the plot. This point is then used as an initial point to feed to fsolve. And yes, fsolve may not succeed if not provided with either ranges or initial points, or it may take much longer, or it may find another solution than the one intended.
If you have no idea in which domain to search for a solution the first could be to use fsolve without initial point or ranges. Otherwise try plotting in various ways.

This variant covers the case where the first set has more than one element and also saves the result instead of just printing it.
It works for e.g.
L:=[{a}, {b, c}, {a, d, f}, {b, d}];
L:=[{a}, {b, c}, { d, f}, {b, d}];
L:=[{a,b}, {b, c}, { a,d, f,b}, {b, d}];

res:='res':
look:=L[1]:
for i from 2 to nops(L) do
 if look subset L[i] then res:=L[i] minus look; break end if
end do;
res;

This variant covers the case where the first set has more than one element and also saves the result instead of just printing it.
It works for e.g.
L:=[{a}, {b, c}, {a, d, f}, {b, d}];
L:=[{a}, {b, c}, { d, f}, {b, d}];
L:=[{a,b}, {b, c}, { a,d, f,b}, {b, d}];

res:='res':
look:=L[1]:
for i from 2 to nops(L) do
 if look subset L[i] then res:=L[i] minus look; break end if
end do;
res;

@OffshoreEngineer Take a look at the following approach, where initially only the conditions at x = 0 are used.

W:=simplify(rhs(dsolve({deqns,inits[1..2]}))) assuming EL>0,IN>0;
indets(W,name) minus indets(deqns,name);
#We now need to impose the boundary conditions at x = 100.
IntegrationTools:-Expand(W);
W1:=collect(%,[_C1,_C3]); #Assuming that the arbitrary constants are indeed labelled 1 and 3
eq1:=eval(W1,x=100)=0;
eq2:=eval(diff(W1,x,x),x=100)=0;
Csol:=solve({eq1,eq2},{_C1,_C3}):
#The final solution is:
wsol:=eval(W1,Csol):
#However, notice that there are lots of integrals
indets(wsol,specfunc(anything,Int));
#Although they are all inert (i.e. Int instead of int) value doesn't help, since Maple cannot evaluate them, which this shows:
wsol-subs(int=Int,value(wsol)); #Result 0

Working with wsol will involve lots of numerical integration and it does not seem to be a good idea at all.

@OffshoreEngineer Take a look at the following approach, where initially only the conditions at x = 0 are used.

W:=simplify(rhs(dsolve({deqns,inits[1..2]}))) assuming EL>0,IN>0;
indets(W,name) minus indets(deqns,name);
#We now need to impose the boundary conditions at x = 100.
IntegrationTools:-Expand(W);
W1:=collect(%,[_C1,_C3]); #Assuming that the arbitrary constants are indeed labelled 1 and 3
eq1:=eval(W1,x=100)=0;
eq2:=eval(diff(W1,x,x),x=100)=0;
Csol:=solve({eq1,eq2},{_C1,_C3}):
#The final solution is:
wsol:=eval(W1,Csol):
#However, notice that there are lots of integrals
indets(wsol,specfunc(anything,Int));
#Although they are all inert (i.e. Int instead of int) value doesn't help, since Maple cannot evaluate them, which this shows:
wsol-subs(int=Int,value(wsol)); #Result 0

Working with wsol will involve lots of numerical integration and it does not seem to be a good idea at all.

@zippon I worked a little more on this.
I seem to be getting correct results all the time if I declare r local:
task := proc (a, i, j) local k,r;
etc.

@zippon I worked a little more on this.
I seem to be getting correct results all the time if I declare r local:
task := proc (a, i, j) local k,r;
etc.

Maybe you could provide examples. It is clearly not always a problem. Example:
A:=[seq](Int(x*exp(-x^(2*n)),x=-infinity..infinity),n=1..10);
value~(A);
evalf(A);

@J4James Here is an ordered version (revised from an earlier version to make it faster):

restart;
Eq1:=diff(u(eta,t),t)=C*diff(u(eta,t),eta$2);
BCs := {u(0,t)=sin(t), u(L,t)=0};
ICs := {u(eta,0)=0};
L:=600:
sol:= C1->pdsolve(eval(Eq1,C=C1),ICs union BCs,numeric,time=t,range=0..L,spacestep=.1);
hC:=.01: ht:=.1: he:=.01: #Increments in C,t, and eta, respectively
Ne:=4: NC:=3: Nt:=2: #Number of increments
A:=Matrix(NC*Nt*Ne,4); #Rows: [C,eta,t,u]
t0:=time():
r:=0:
for i from 1 to NC do
   q:=sol(i*hC);
   for j from 1 to Nt do
     v:=q:-value(t=j*ht);
     for k from 1 to Ne do
       r:=r+1;
       A[r,..]:=subs(v(he*k),Vector[row]([i*hC,eta,t,u(eta,t)]))
     end do
   end do
end do;
time()-t0;
#interface(rtablesize=infinity);
A;

@J4James Here is an ordered version (revised from an earlier version to make it faster):

restart;
Eq1:=diff(u(eta,t),t)=C*diff(u(eta,t),eta$2);
BCs := {u(0,t)=sin(t), u(L,t)=0};
ICs := {u(eta,0)=0};
L:=600:
sol:= C1->pdsolve(eval(Eq1,C=C1),ICs union BCs,numeric,time=t,range=0..L,spacestep=.1);
hC:=.01: ht:=.1: he:=.01: #Increments in C,t, and eta, respectively
Ne:=4: NC:=3: Nt:=2: #Number of increments
A:=Matrix(NC*Nt*Ne,4); #Rows: [C,eta,t,u]
t0:=time():
r:=0:
for i from 1 to NC do
   q:=sol(i*hC);
   for j from 1 to Nt do
     v:=q:-value(t=j*ht);
     for k from 1 to Ne do
       r:=r+1;
       A[r,..]:=subs(v(he*k),Vector[row]([i*hC,eta,t,u(eta,t)]))
     end do
   end do
end do;
time()-t0;
#interface(rtablesize=infinity);
A;

@geri23 It may be worth your while to formulate your problem again from the very start and then post that as a new question in MaplePrimes. That way other people get a chance to participate in the discussion. The present thread is likely read by only you and me.