@digerdiga Mistakingly in indicialeq I took the second argument to mean the name of the variable to be used in the indicial equation. Since I wanted the name r, I wrote indicialeq(ode2, r, 0, X(y));
And that produces r^3-3*r^2+2*r = 0.
However, the second argument has to be the name of the independent variable in the ode, in this case y.
My mistake.
In my comment above I deleted the passages you are quoting.

############################
The roots of the indicial equation are 1, 5/2, and 4. Two of them differ by an integer.
That is well known to complicate matters (see any book that covers series solutions, e.g. Boyce & diPrima or Coddington & Levinson). But the usual procedure works for 5/2 and for the larger of the two roots (i.e. 4).
For 5/2 we can proceed as I did in a comment to your original question:
eval(ode2,X(y)=y^(5/2)*u(y));
ode21:=collect(expand(%/y^(5/2)),diff,factor);
f1:=diffeqtorec(ode21, u(y), a(n));
sol1:=rsolve({f1,a(0)=a0,a(1)=a1},a(n),makeproc);
X1:=y^(5/2)*add(sol1(i)*y^i,i=0..3); #Truncated
#For the root 4 you can do the same thing with 4 taking the place of 5/2.
#For the root 1 which differs from 4 by an integer, there will be problems like the ones discussed earlier about division by zero. However, as it happens, there is the polynomial solution we found earlier.

@digerdiga The command
dsolve(ode2, X(y), 'formal_series', y = 0);
ignores the last argument since it allows considers 'coeffs'=whatever only.
Notice that you get expansions in powers of y, (y+1), and (y-1/z), in accordance with the coefficient of X'''(y), which is y^3*(y-1)*(y*z-1).

@digerdiga The command
dsolve(ode2, X(y), 'formal_series', y = 0);
ignores the last argument since it allows considers 'coeffs'=whatever only.
Notice that you get expansions in powers of y, (y+1), and (y-1/z), in accordance with the coefficient of X'''(y), which is y^3*(y-1)*(y*z-1).

@digerdiga The algorithm solves for a(n+2) (try showstat(sol1);)
So
f1:=op(solve(f[1],{a(n+2)}));
#You see that necessarily you get a numeric exception when doing:
for nn from 0 to 2 do eval(f1,n=nn) end do;

@digerdiga The algorithm solves for a(n+2) (try showstat(sol1);)
So
f1:=op(solve(f[1],{a(n+2)}));
#You see that necessarily you get a numeric exception when doing:
for nn from 0 to 2 do eval(f1,n=nn) end do;

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

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