@digerdiga There is a polynomial solution:

dsolve(ode2,X(y),'formal_solution');

#That can also be found by:
sol1:=rsolve({f[1],a(0)=a0,a(1)=a1},a(n),makeproc);
X(y)=add(sol1(i)*y^i,i=0..3);
odetest(%,ode2);
#You see that a0 must be 0.
#To see what is going on you could try to determine a(
for nn from 0 to 4 do eval(f[1]=0,n=nn) end do;
#Non-polynomial solutions:
sol2:=rsolve({f[1],a(4)=1,a(5)=0},a(n),makeproc);
X(y)=add(sol2(i)*y^i,i=4..8); ##Truncated of course
sol3:=rsolve({f[1],a(4)=0,a(5)=1},a(n),makeproc);
X(y)=add(sol3(i)*y^i,i=4..8);

What is the question?

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

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