As far as dsolve is concerned your syntax is not quite right. Try this:

restart;
mystery := proc (z) evalf((1+I)*z) end proc;
F1 := proc (N, t, Y, YP) local tmp;
   tmp:=mystery(Y[1]+I*Y[2]);
   YP[1] := Re(tmp);
   YP[2] := -Im(tmp)
end proc;
curve := dsolve(numeric,procedure = F1, number = 2, procvars = [x(t), y(t)], start = 0, initial = Array([-.2,.3]));
plots:-odeplot(curve,[x(t),y(t)],-1..1);
curve(1.2345);
#The last output with the .0*I terms hints at the possible problem DEplot is having.
#At the moment I only have access to Maple 12, but there it seems that the parameters option doesn't work when the parameters don't appear in the system, in this case in F1.

dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-x(t),x(0)=1,y(0)=0};

dsn := dsolve(dsys,numeric,output=listprocedure);

X,Y:=op(subs(dsn,[x(t),y(t)]));
evalf(Int(X(t),t=0..1));   
#For your simple system you could have found the integral from the second equation much simpler:
Y(0)-Y(1);

If the system is more complicated you can still use the latter method simply by introducing z(t) by
diff(z(t),t) = x(t) with z(0) = 0. Then z(1) = int(x(t),t=0..1):

restart;
dsys2 := {diff(x(t),t)=y(t),diff(y(t),t)=-sin(x(t)),diff(z(t),t)=x(t),x(0)=1,y(0)=0,z(0)=0};
dsn2 := dsolve(dsys2,numeric,output=listprocedure);
X,Y,Z:=op(subs(dsn2,[x(t),y(t),z(t)]));
Z(1);
#Compare with:
dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-sin(x(t)),x(0)=1,y(0)=0};
dsn := dsolve(dsys,numeric,output=listprocedure);
X,Y:=op(subs(dsn,[x(t),y(t)]));
evalf(Int(X(t),t=0..1));   

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