If your equation really has the form
ode:=diff(y(x),x)=w(y(x))
for some known function w, then with Y:=eval(y(x),dsol8) you could simply do
Y1:=w@Y

You could use a loop instead of sequence:

restart;
ode:=diff(x(t),t)=x(t)^2;
dsolve({ode,x(0)=1},numeric);
rhs(res(.5)[2]);
rhs(res(1.1)[2]); #error
seq(rhs(res(m/10)[2]),m=0..40); #error
#Use a loop and break after the first error. No point in trying higher m-values.
S:=NULL:
for m from 0 to 40 do
   try
     S:=S,rhs(res(m/10)[2])
   catch:
     break
   end try
end do;
 
S;
m;
#If you don't mind the error then here is the short version:
S:=NULL:
for m from 0 to 40 do S:=S,rhs(res(m/10)[2]) end do:
S;
m;

As another idea you could do:

resA:=dsolve({ode,x(0)=1},numeric,output=Array([seq(m/10,m=0..40)]));
A:=resA[2,1];
plot(A);
A[11,2];
whattype(%);

Finally a method that allows you to use seq (using res from above):

X:=proc(t) local r; try r:=rhs(res(t)[2]) catch: r:=NULL end try; r end proc;
X(.5);
X(1);
seq(X(m/10),m=0..40);

There are several issues, but take a look at this modified example:

restart;
a := 3;
U := u(r, t);
wave := a^2*(diff(U, r$2)+diff(U, r)/r) = diff(U,t$2);
bcs := u(1, t) = 0, u(0,t)=1; #added a boundary condition
ics := D[2](u)(r,0)=piecewise(r<1/2,-2,0) ,u(r, 0) = 0; #Rewrote initial conditions
s := pdsolve(wave, {bcs, ics},numeric,time=t,range=0..1); #Solving numerically
s:-animate(t=5);

If I understand your worksheet correctly, the example mentioned on the help page for events may be what you need.
?dsolve,numeric,events
Here I choose a simple equation with two discrete variables u(t) and n(t).
restart;
eq:=diff(y(t),t,t)+u(t)+y(t)=sin(n(t)*t);
res:=dsolve({eq,y(0)=0,D(y)(0)=1,u(0)=0,n(0)=1},numeric,output=listprocedure,
  discrete_variables=[u(t)::float,n(t)::float],
  events=[[t=[0,1/2],u(t)=y(t)],[t=[0,1/2],n(t)=2*t]]);
plots:-odeplot(res,[t,y(t)],0..30);
plots:-odeplot(res,[t,n(t)],0..30);
plots:-odeplot(RES,[t,u(t)],0..30);

restart;
plot3d([sqrt(x^2+y^2),sqrt(1-x^2-y^2)],x=-1..1,y=-1..1);
#Somewhat rough, so for plotting purposes we change to cylindrical coordinates
?plot3d,coords
addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);
plot3d([r,sqrt(1-r^2)],r=0..1,theta=0..2*Pi,coords=z_cylindrical,scaling=constrained,transparency=[0,.4]);
#Using spherical coordinates r, theta, phi in the computation:
Int(Int(Int(exp(r^3)*r^2*sin(theta),r=0..1),theta=Pi/4..Pi/2),phi=0..2*Pi);
value(%);

Int(Int(sqrt(x^2+y^2),x=y..sqrt(2*y-y^2)),y=0..1);
value(%);
simplify(%);
#In polar coordinates:
Int(Int(r^2,r=0..2*sin(theta)),theta=0..Pi/4);
value(%);

Before assigning to m you can try isolating the highest derivatives:
solve({Eq1,Eq2},{diff(f(eta),eta,eta),diff(theta(eta),eta,eta)});
simplify(%) assuming diff(f(eta),eta)<0; #or f' > 0 if you like

Now you got an equation of the form f'' = expr, where expr contains f' raised to a power less than 1. That is likely to cause problems if f' becomes zero in the interval, which unfortunately it must, since you require f(-1) = 0 and f(1) = 0.
Just think of the nonuniqueness properties of y' = y^(1/2), y(0)=0.

Notice that the integral over z has to be restricted to z = 0..1.
restart;
with(Statistics):
dummy := sqrt((1+y*e)/(y*e));
dummy1 := sqrt(1/(1+y*e));
Q := RandomVariable(Normal(0, dummy));
R := RandomVariable(Normal(y*e*z/(1+y*e), dummy1));
y := 1;
e := 1/5;
k := 1;
l:=1/2;
a:= 5/4;
combine(Int(k*r*PDF(R,r),r=-infinity..x),power);
J:=value(%);
limit(J,x=-infinity);
factor(diff(J,x));
#So J < 0 for all x < 0 and all z.
#Furthermore J is an increasing function of x  for x > 0. For z = 0 J < 0 for all x.
#For z > 0 there is a solution since
limit(J,x=infinity);
#This means that you have to restrict the integral over z from -1..1 to 0..1.
#An illustration:
plots:-animate(plot,[J,x=-5..5],z=-1..1);
#Search interval can be restricted to x=0..infinity since J < 0 for all x < 0:
F:=unapply('fsolve'(J=0,x=0..infinity),z,numeric);
F(1);
F(0); #Output NULL
F(.5);
F(0.0000001);
evalf(Int(F(z)*PDF(Q,z),z=0..1));
#Result 0.2829758111 (Digits at 10)

With symbolic input I'm sure you don't really want to see the eigenvalues.
A:=Matrix([[-X,0,-beta2*S,alpha],[beta*IB,-K,0,0],[beta2*IM,epsilon,Y,0],[0,gamma,sigma,-W]]);
#You get the eigenvalues by
Lambda:=LinearAlgebra:-Eigenvalues(A):
#However, the result (which I didn't wait for) takes up a lot of space.
#Consider
LinearAlgebra:-CharacteristicPolynomial(A,lambda);
#This polynomial is not really any simpler than the general one:
p4:=lambda^4+add(a[i]*lambda^i,i=0..3);
solve(p4=0,lambda);
allvalues(%);
length(%);
#Imagine substituting for a[i] the coefficients of the characteristic polynomial!

Leaving the ode somewhat abstract can help:

restart;
Eq2 := diff(z(x), x, x) = p(x)*f(diff(z(x), x));
p:=x->piecewise(x<10,2,x<11,10,2);
res:=dsolve({Eq2,z(0)=0,D(z)(0)=0});
eval(res,f=(u->sqrt(1+u^2))); #Now making f concrete
map(value,rhs(%));
evalindets(%,RootOf,allvalues);
simplify(value(%));

I'm not very satisfied with this  solution, but present it anyway:
restart;
poly:=2*x^3-4*x^2-22*x+24;
solve(poly=0,x); #So we know the results, which also can be written (apparently):
r1:=2/3 * 37^(1/2) * sin(1/6 * Pi+1/3 * arccos(55/1369 * 37^(1/2)))+2/3;
r2:=-1/3 * 37^(1/2) * (sin(1/6 * Pi+1/3 * arccos(55/1369 * 37^(1/2)))+3^(1/2)*sin(1/3 * Pi-1/3 * arccos(55/1369 * 37^(1/2))))+2/3;
r3:=-1/3 * 37^(1/2) * (sin(1/6 * Pi+1/3 * arccos(55/1369 * 37^(1/2)))-3^(1/2) * sin(1/3 * Pi-1/3 * arccos(55/1369 * 37^(1/2))))+2/3;
#Easily checked that they are correct:
collect((x-s1)*(x-s2)*(x-s3),x,expand);
combine(expand(r1*r2*r3));
combine(expand(r1+r2+r3));
combine(expand(r1*r2+r1*r3+r2*r3));
#Now more directly:
u:=arccos(55/1369 * 37^(1/2));
subs(v=w/3,cos(3*v)=expand(cos(3*v)));
Rc:=solve(%,{cos(w/3)});
subs(v=w/3,sin(3*v)=expand(sin(3*v)));
Rs:=subs(Rc[1],solve(%,{sin(w/3)}));
Rc1,Rs:=op(simplify(subs(w=u,[Rc[1],Rs])));
subs(Rc1,Rs,expand(r1));
subs(Rc1,Rs,expand(r2));
subs(Rc1,Rs,expand(r3));

I take here an example from the paper by Amendola et al. (http://arxiv.org/abs/gr-qc/0612180), which you gave as a reference.

On page 6 they mention as an example f(R) = R+alpha*R^(-n), which they mention corresponds to m(r) = -n*(1+r)/r with n a constant. This example is mentioned again on page 7 under the discussion of P2, this time though with f(R) = R+alpha*R^n, which similarly corresponds to m(r) = n*(1+r)/r.
restart;
f := (x, y, z)-> -1-z-3*y+x^2-x*z;
g := (x, y, z)-> x*z/m(z/y)-2*z*y+4*y+x*y;
h := (x, y, z)-> -x*z/m(z/y)-2*z^2+4*z;
m:=s->n*(1+s)/s; #page 7 under P2 = (-1,0,0)
res := solve({f(x, y, z) = 0, g(x, y, z) = 0, h(x, y, z) = 0}, {x, y, z});
pts:=map(subs,[res],[x,y,z]);
Looking at (-1,0,0), which wasn't found by Maple
f(-1,0,0); #OK zero
map(normal,g(-1,y,z)); #This has no limit as (y,z) -> (0,0)
map(normal,h(-1,y,z)); #This doesn't either
#The problem is that z^2/(y+z) doesn't have a limit if y and z are allowed different signs:
#Approach (y,z) = (0,0) along y = - z + z^3 as an example.
#But if y and z must have the same sign (positive, say) then z^2/(z+y) <= (z/(z+y))*z <= z -> 0 as (y,z) -> (0,0) through positive values. However, since the first critical point P1 is (0,-1,2) it appears that different signs of y and z are acceptable. Thus P2 is not a critical point for this system.

#No problem with the other points:
J := unapply(VectorCalculus:-Jacobian([f(x, y, z), g(x, y, z), h(x, y, z)], [x, y, z]), x, y, z):
use ev=LinearAlgebra:-Eigenvalues in
  for p in pts do
  ev(J(op(p)))
  end do
end use;

Maybe converting to LegendreP could help:

restart;
chi:=Pi/2-theta:
Vsi:=-(1/sqrt( 1+epsilon - cos(chi)*cos(phi-Pi/3)) + 1/sqrt( 1 +epsilon- cos(chi)*cos(phi-Pi))+ 1/sqrt(1+epsilon-cos(chi)*cos(phi+Pi/3)));
A:=Int(Vsi*SP, [theta=0..Pi,phi=0..2*Pi]);
SP1:=SphericalY(3,2,theta,phi)*conjugate(SphericalY(2,1,theta,phi))*sin(theta);
convert(SP1,LegendreP);
SP2:=simplify(%)  assuming real;
#This can be written as follows:
SP3:=-(15/16*I)*exp(I*phi)*cos(theta)^2*sqrt(7)*sin(theta)^4/Pi;
#Just a check of orthogonality:
int(SP3, [theta=0..Pi,phi=0..2*Pi]);
#It may or may not speed up calculations to split in real and imaginary parts:
SPr,SPi:=op(evalc([Re,Im](SP3)));
evalf(subs(SP=SPr,epsilon=0.001,A));
evalf(subs(SP=SPi,epsilon=0.001,A));

Although RootOf's may pop up, your main problem seems to be that you are missing a multiplication sign in

0.5*pi(Rmax^4-x^4)

This should be
0.5*Pi*(Rmax^4-x^4)

Notice that I also changed pi to Pi.

If after correcting this you still get a RootOf then try allvalues.
Finally: It is much easier to get a good response if you post the code as text or as an uploaded worksheet. That way people are more likely to try your code.

The definition in the help page for FourierTransform gives the transform Z of z as

eq:=Z[i] = sqrt(1/N)*Sum(z[j]*exp(-(2*I)*Pi*(i-1)*(j-1)/N), j = 1 .. N);# i = 1 .. N
#so
eval(eq,i=1);
#which implies for your v that Z[1] is real.
#Further for N even (as in your case)
eval(eq,i=N/2+1);
#which implies that Z[N/2+1] can be written as seen here explicitly when N=6: 
value(eval(%,N=6));
#This means that in your case Z[4] will be purely imaginary.
#Doing the computation in your case:
v:=<1+2*I,5+3*I,9+4*I,9-4*I,5-3*I,1-2*I>;
N:=6:
Z := Vector(N,i->evalf(sqrt(1/N)*add(v[j]*exp(-(2*I)*Pi*(i-1)*(j-1)/N), j = 1 .. N)));
#which agrees with
DiscreteTransforms:-FourierTransform(v);
####
However, if you take a look at
http://en.wikipedia.org/wiki/Fast_Fourier_transform
you see mentioned that if the input is a real vector then the transform satisfies:
conjugate(Z[i]) =Z[N-i+2]; #i=2..N
Now try
w:=<5,1+2*I,5+3*I,9,5-3*I,1-2*I>;
and
DiscreteTransforms:-FourierTransform(w);
DiscreteTransforms:-InverseFourierTransform(w);