@Markiyan Hirnyk So you are back to your old self.

@Markiyan Hirnyk Surely it is sufficient that M be negative definite for the system x' = M.x to be (asymptotically) stable. But negative (and positive) definiteness implies in particular that all the eigenvalues are real. That it is surely not necessary for the system to be stable.

@Markiyan Hirnyk Surely it is sufficient that M be negative definite for the system x' = M.x to be (asymptotically) stable. But negative (and positive) definiteness implies in particular that all the eigenvalues are real. That it is surely not necessary for the system to be stable.

Why do you want an answer to this (besides general curiosity)?
If the question has to do with stability of the system x' = M.x then a more relevant question is: What are the necessary and sufficient conditions on m and q for all eigenvalues of M to have negative real part?
For that question you can apply the Routh-Hurwitz criterion on the characteristic polynomial of M.
See ?PolynomialTools:-Hurwitz

Why do you want an answer to this (besides general curiosity)?
If the question has to do with stability of the system x' = M.x then a more relevant question is: What are the necessary and sufficient conditions on m and q for all eigenvalues of M to have negative real part?
For that question you can apply the Routh-Hurwitz criterion on the characteristic polynomial of M.
See ?PolynomialTools:-Hurwitz

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